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Abstract 

The  United  States  Air  Force  and  the  Department  of  Defense  is  increasingly  inter¬ 
ested  in  nanomaterials.  To  stndy  these  materials,  one  needs  to  measure  the  mechanics 
of  materials  on  the  nanoscale.  Over  the  past  few  decades  the  atomic  force  microscope 
(AFM)  has  been  used  in  various  methods  to  establish  local  surface  properties  at  the 
nanoscale.  In  particular,  surface  elasticity  measurements  are  crucial  to  understanding 
nanoscale  surface  properties.  Problems  arise,  however,  when  measuring  soft  surfaces 
such  as  polymers  and  biological  specimens,  because  these  materials  have  a  more  com¬ 
plex  viscoelastic  response. 

This  research  focuses  on  modeling  an  AFM  dynamic  nanoindentation  experiment 
intended  to  characterize  near-surface  viscoelastic  material  parameters.  The  exper¬ 
iment  uses  an  AFM  in  dynamic  contact  mode  with  a  polymer  surface  to  gather 
frequency  dependent  amplitude  and  phase  data.  A  three-dimensional,  dynamic  vis¬ 
coelastic  model  of  the  AFM  and  surface  interaction  is  developed  and  then  analytically 
solved  in  the  linear  approximation  under  appropriate  physical  assumptions  based  on 
the  physics  of  the  AFM  experimental  setup.  As  an  illustrative  application,  the  an¬ 
alytical  solution  is  coupled  with  experimental  data  from  a  polystyrene  material  to 
ascertain  surface  material  properties  at  the  nanoscale.  Our  solution  allows  the  direct 
calculation  of  the  storage  and  loss  modulus  from  experimental  data. 
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AN  ANALYTICAL  MODEL  OF  NANOMETER  SCALE  VISCOELASTIC 
PROPERTIES  OF  POLYMER  SURFACES  MEASURED  USING  AN  ATOMIC 

FORCE  MICROSCOPE 

I.  Introduction 

The  United  States  Air  Force  and  the  Department  of  Defense  is  increasingly  inter¬ 
ested  in  nanomaterials  [17].  In  order  to  study  these  materials,  one  needs  to  measure 
the  mechanics  of  materials  on  the  nanometer  scale.  With  technology’s  increasing  fo¬ 
cus  on  the  atomic  and  molecular  level,  the  importance  of  nanoscale  surface  physics 
and  chemistry  cannot  be  understated.  Since  the  invention  of  the  atomic  force  micro¬ 
scope  (AFM)  in  1986  [9],  there  has  been  substantial  growth  in  its  applications.  In 
particular,  the  AFM  has  been  used  to  establish  local  surface  properties  such  as  elastic 
modulus,  adhesion,  surface  friction,  viscosity,  hardness,  energy  dissipation,  and  glass 
transition  temperature  [3,  10,  19,  21,  23,  40,  45,  90,  49,  59,  68,  76,  86,  87,  88,  91,  99]. 
Determining  nanoscale  surface  properties  is  a  challenge  at  the  atomic  level  because 
surface  materials  behave  different  from  the  bulk  material.  Forces  such  as  adhesion  and 
van  der  Waals  bond  interactions,  which  are  usually  unsubstantial  at  the  macro  scale, 
now  become  dominant  influences.  In  addition,  problems  arise  when  dealing  with 
softer  surfaces  such  as  polymers  and  biological  specimens.  One  reason  is  because 
these  materials  have  more  complex  viscoelastic  and  even  plastic  responses.  Another 
important  reason  is  that  these  materials  can  have  anisotropic  responses. 

This  research  focuses  on  modeling  an  AFM  dynamic  nanoindentation  experiment 
on  a  viscoelastic  surface.  In  the  experiment,  an  AFM  in  dynamic  contact  mode  with 
a  polymer  surface  is  used  to  gather  frequency  dependent  amplitude  and  phase  data. 
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A  mathematical  model  of  the  experiment  is  created  and,  under  several  appropriate 
physical  assumptions,  an  analytical  solution  in  the  linear  approximation  is  achieved 
for  the  penetration  depth  of  the  AFM  tip  into  the  material.  Another  model  is  created 
that  relates  the  raw  experimental  data  to  the  depth  of  penetration  of  the  AFM  tip. 
These  two  models  are  then  equated  to  ascertain  an  analytical  solution  allowing  the 
calculation  of  frequency  dependent,  viscoelastic  material  properties  from  the  experi¬ 
mental  data.  As  an  illustration,  viscoelastic  properties  for  a  polystyrene  material  are 
calculated  and  presented  using  the  analytical  solution. 

This  chapter  is  divided  into  three  sections.  The  hrst  section  contains  the  motiva¬ 
tion  for  this  research  along  with  other  research  that  has  been  performed.  The  second 
section  gives  a  brief  overview  of  the  experiment  that  is  to  be  modeled.  Finally,  we 
detail  what  will  be  presented  in  the  remaining  chapters  of  this  document. 

1.1  Motivation 

The  interest  of  the  United  States  Air  Force  and  the  Department  of  Defense  in 
nanomaterials  is  of  key  importance  in  many  areas.  In  particular,  high-altitude  long- 
endurance  ISR  airships,  prompt  theater-range  ISR/strike  systems,  direct  forward  air 
delivery  and  resupply,  energy-efficient  partially  buoyant  cargo  airlifters,  fuel-efficient 
hybrid  wing-body  aircraft,  and  hyperprecision  low-collateral  damage  munitions  [17]. 
In  order  to  help  further  research  in  these  areas,  nanoscale  viscoelastic  material  prop¬ 
erties  are  of  particular  importance.  These  properties  are  calculated  using  both  math¬ 
ematical  theory  and  experiment.  Though  the  experiments  widely  vary,  their  mathe¬ 
matical  approach  remains  consistent,  that  is  the  movement  of  the  indentation/loading 
system,  whether  an  AFM  or  a  nanoindenter,  is  often  modeled  as  a  spring-mass  sys¬ 
tem  [5,  97,  98].  In  addition,  the  contact  model  of  the  indentation/loading  system  with 
the  material  surface  is  the  Sneddon  solution  [84]  extended  to  the  viscoelastic  case  using 
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the  method  of  functional  equations  proposed  by  Lee  and  Radok  [52],  Sneddon’s  solu¬ 
tion  is  derived  from  a  rigid  axisymmetric  indenter  that  is  indented  into  a  half-space 
composed  of  a  homogeneous,  isotropic,  linearly  elastic  material.  Sneddon’s  solution 
is  based  on  the  quasi-static  equilibrium  equations  and  provides  the  penetration  depth 
in  terms  of  the  indenter’s  geometry  as  well  as  the  load  imposed  on  the  indenter  in 
terms  of  its  geometry.  Using  the  method  of  functional  equations,  Sneddon’s  solu¬ 
tion  is  extended  to  the  linearly  viscoelastic  case  provided  that  the  contact  area  of 
the  indenter  increases  monotonically  with  time.  This  extended  solution  provides  a 
mathematical  way  of  relating  monotonically  increasing  and/or  constant  load  histo¬ 
ries  to  the  complex  stiffness  of  a  material  given  a  chosen  viscoelastic  model  [67,  81]. 
The  extended  solution,  along  with  a  particular  viscoelastic  model,  can  also  be  used 
in  specihc  testing  procedures  such  as  the  load  relaxation  test  or  the  creep  test,  to 
produce  the  relaxation  modulus  or  creep  compliance,  respectively  [15,  96,  97]. 

We  model  the  dynamic  nanoindentation  experiment  in  two  pieces  which,  are  cou¬ 
pled  together  using  the  appropriate  boundary  conditions.  The  two  pieces  of  the 
experimental  setup  are  the  AFM  and  the  material.  The  AFM  is  modeled  as  a  one¬ 
dimensional  spring-mass  system  with  three  driving  forces.  The  driving  forces  consist 
of  the  static  and  sinusoidal  forcing  provided  by  the  AFM  system,  and  the  force  on 
the  tip  caused  from  the  stress  within  the  surface.  The  latter  force  serves  to  couple 
the  AFM  with  the  material.  The  modeling  of  the  AFM  as  a  spring-mass  system  is 
similar  to  current  AFM  modeling  theory  [53,  95]. 

The  material  is  modeled  based  on  continuum  mechanics  using  a  three-dimensional, 
linear  viscoelastic  material  model.  The  generated  equations  are  the  elastodynamic 
equilibrium  equations,  i.e.,  they  contain  an  inertial  term.  This  clearly  differs  from 
the  current  theories  which  are  based  upon  the  quasi-static  equilibrium  equations  [15, 
56,  58,  67,  81,  96].  An  axisymmetric  indenter,  representing  the  AFM  tip,  is  applied 
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to  the  surface  of  the  material.  The  application  of  the  indenter  is  such  that  there  is  a 
static  indentation  superposed  with  a  much  smaller  sinusoidally  varying  indentation. 
This  represents  the  nature  of  the  experiment  in  that  the  driving  force  from  the  system 
contains  a  static  loading  which  is  far  greater  than  the  dynamic  loading.  The  related 
equations  are  then  analytically  solved  in  the  linear  approximation  under  the  assump¬ 
tion  that  the  static  loading  far  exceeds  the  dynamic  loading.  The  analytic  solution  is 
in  terms  of  the  penetration  depth  of  the  indenter  into  the  material. 

A  model  that  relates  the  raw  experimental  data  to  the  penetration  depth  of  the 
AFM  tip  is  then  equated  to  the  analytical  solution  from  the  material  model.  This 
allows  for  the  calculation  of  the  frequency  dependent  storage  and  loss  modulus  based 
on  a  generic  viscoelastic  polynomial  model.  Our  solution  allows  the  direct  calcula¬ 
tion  of  the  storage  and  loss  modulus  from  experimental  data.  A  specihc  viscoelastic 
model  can  then  be  chosen  to  £t  the  data  as  accurately  as  possible  in  order  to  attain 
specihc  material  parameters.  This  differs  from  current  solutions  in  that  current  so¬ 
lutions  have  to  choose  a  viscoelastic  model  before  the  storage  and  loss  modulus  are 
calculated  [15,  67,  81,  96].  Thus,  our  model  allows  for  a  more  varied  approach  to 
choosing  a  viscoelastic  material  model. 

This  concludes  the  overview  of  the  motivation  for  our  research.  The  next  section 
gives  an  overview  of  the  experimental  setup  and  procedure. 

1.2  Experimental  Setup 

The  experiment  we  are  concerned  with  uses  an  AFM  as  a  dynamic  indentor  in 
order  to  glean  nanomechanical  properties  of  a  sample  surface.  The  representation 
of  the  experimental  setup  can  be  seen  in  Figure  1.  In  the  experiment,  the  AFM 
cantilever  tip  is  moved  toward  the  sample  surface  by  expansion  of  the  piezoelectric 
crystal  (Piezo).  The  considered  sample  consists  of  a  polystyrene  layer  bonded  to  a 
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Error  L-l-A 


Figure  1.  A  simplified  model  of  the  experimental  setup. 

silicon  base.  The  silicon  base  rests  atop  a  piezo.  A  voltage  is  applied  to  the  piezo, 
which  causes  the  crystal  to  expand  thus  raising  the  sample’s  surface  into  the  AFM 
tip.  This  voltage  is  the  DC  component  of  the  piezo  drive  signal.  Since  the  surface 
is  viscoelastic,  there  is  some  initial  creep,  but  after  a  set  amount  of  time,  the  system 
approaches  an  equilibrium  position. 

A  laser  beam  is  continuously  reflected  off  the  tip  end  of  the  cantilever  and  onto  a 
detector.  The  initial  voltage  measured  at  the  detector  is  called  the  preload  voltage,  or 
the  static  load.  The  preload  voltage  is  considered  the  setpoint  of  the  experiment.  It 
is  the  voltage  which  corresponds  to  a  cantilever  position  around  which  the  cantilever 
will  be  oscillating.  This  voltage  is  considered  to  be  the  initial  start  point  or  zero 
point. 

A  manual  bias  is  set  on  the  detector.  The  voltage  of  this  manual  bias  is  called 
the  modulation  voltage  or  forcing  amplitude,  and  its  phase  is  the  forcing  frequency  or 
tapping  frequency.  This  forcing  never  changes  and  is  also  our  reference  amplitude  and 
phase.  A  closed-loop,  computer  controlled  feedback  system  modihes  the  modulation 
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voltage  and  phase  in  an  attempt  to  maintain  the  setpoint  within  the  modulation 
voltage  bias.  The  total  signal  that  is  input  to  the  piezo  is  called  the  piezo  drive 
signal.  A  lock-in-amplifier  (LIA),  which  we  call  the  AC  component  of  the  piezo  drive 
signal,  measures  the  exact  piezo  voltage  input  and  phase.  The  piezo  LIA  records  a 
continuous  measurement  for  all  tapping  frequencies. 

The  piezo  drive  signal  forces  the  piezo  to  oscillate  at  a  fixed  frequency  and  am¬ 
plitude,  causing  the  movement  of  the  surface  relative  to  its  initial  position  to  be  of 
an  oscillatory  nature.  Since  the  piezo  is  a  crystal,  it  takes  a  set  amount  of  time  to 
contract  or  expand  to  the  voltage  that  is  applied  to  it.  This  delay  in  contraction  or 
expansion  is  called  the  piezo  delay  time,  which  is  in  this  case  frequency  dependent. 
Consequently,  the  oscillatory  nature  of  the  surface’s  relative  position  occurs  at  the 
piezo  drive  signal,  but  it  is  phase  delayed  by  the  piezo  delay  time. 

As  the  surface’s  relative  position  oscillates,  it  interacts  with  the  AFM  tip  and 
causes  the  tip  to  move  in  an  oscillating  pattern.  This  tip  movement  is  not  at  the 
same  amplitude  or  phase  as  the  piezo  drive  signal  nor  is  it  at  the  same  amplitude  and 
phase  as  the  modulation  voltage.  The  interaction  between  the  soft  sample  surface  and 
the  harder  AFM  tip  creates  a  smaller  voltage  amplitude  and  an  additional  delay  time, 
called  the  viscoelastic  retardation  time,  which  is  manifested  in  the  AFM  cantilever 
tip’s  relative  motion.  The  difference  between  the  actual  AFM  tip  movement  and  the 
piezo  modulation  voltage  is  measured  by  a  lock-in-amplifier.  We  call  this  the  error 
LIA.  The  error  LIA  continuously  measures  the  AC  amplitude  and  phase  difference 
for  all  tapping  frequency. 

This  procedure  is  performed  through  a  sweep  of  tapping  frequencies  from  1000 
Hz  to  1  Hz,  allowing  equilibrium  to  be  reached  at  each  frequency.  After  sweeping 
through  all  the  frequencies,  the  process  is  started  over  with  a  new  preload  voltage. 
Between  ten  and  fifteen  different  preload  voltages  are  used,  ranging  from  0.25  V  to 


6 


5.4  V.  The  same  experiment  is  also  performed  using  polystyrene  samples  of  differing 
thickness:  30  nm,  70  nm,  and  220  nm. 

This  concludes  the  discussion  on  the  experimental  setup.  Our  focus  will  be  to 
mathematically  model  the  above  procedure  as  closely  as  possible.  Through  our  pro¬ 
ceedings  we  will  make  simplifying  assumptions  based  on  the  physics  of  the  experiment 
in  order  to  achieve  an  analytical  solution  but,  we  will  maintain  the  pertinent  physics 
needed  to  fully  describe  the  experiment.  The  next  section  of  this  chapter  will  detail 
the  contents  contained  in  the  remainder  of  this  dissertation. 

1.3  Dissertation  Contents 

Chapter  one  contains  an  introduction  to  nanoscale  material  modeling,  as  well  as 
a  description  of  the  experiment  and  the  mathematical  model  that  we  will  produce. 
Chapter  two  includes  a  description  of  AFM  modeling  theories  and  the  formulation  of 
a  one-dimensional  spring-mass  model  that  represents  the  AFM  cantilever.  In  Chapter 
three  we  present  material  modeling  theories  and  create  a  three-dimensional,  viscoelas¬ 
tic  material  model.  The  one- dimensional  AFM  model  and  the  three-dimensional  ma¬ 
terial  model  are  then  nondimensionalized.  In  Chapter  four,  we  present  necessary 
background  material  and  then  proceed  to  make  physical  assumptions  that  allow  us 
to  simplify  the  three-dimensional,  viscoelastic  material  model,  as  well  as  the  one¬ 
dimensional  AFM  model.  This  reduced  three-dimensional  system  is  then  analytically 
solved  in  the  linear  approximation,  and  analyzed.  Chapter  hve  includes  the  creation 
of  a  model  that  describes  how  the  data  was  collected.  This  new  model  is  then  coupled 
to  the  reduced  three-dimensional  model;  and  in  conjunction  with  the  experimental 
data,  a  method  for  determining  viscoelastic  material  properties  is  produced.  Finally, 
Chapter  six  draws  conclusions  and  investigates  areas  for  future  work. 
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II.  The  AFM 


Atomic  force  microscopy  techniques  have  been  around  since  its  invention  in  1986  [9]. 
Various  techniques  have  been  developed  ranging  from  measuring  surface  details  to 
modifying  surface  structures.  This  chapter  focuses  on  aspects  of  AFM  modeling  the¬ 
ory  and  how  the  AFM  can  be  coupled  to  the  surface  it  is  in  contact  with.  The  chapter 
ends  by  choosing  an  AFM  model  with  the  appropriate  surface  contact  model. 

2.1  AFM  Modeling  Theories 

The  process  of  modeling  how  dynamic  atomic  force  microscopy  techniques  are 
used  to  determine  surface  properties,  allows  us  to  gain  a  better  understanding  of  the 
mechanical  properties  of  the  surface  at  the  nanoscale.  An  AFM  is  an  instrument 
used  to  produce  high  resolution,  three-dimensional  images  of  sample  surfaces  on  a 
nanometer  scale.  Using  a  scanning  technique  the  AFM  is  capable  of  measuring  small 
forces  between  the  AFM  tip  and  the  sample  surface. 

The  AFM  tip  is  mounted  on  a  cantilevered  beam  having  a  very  small  mass.  As  the 
tip  is  scanned  across  the  surface,  the  cantilevered  beam  flexes  based  on  the  surface 
topology.  The  flexing  motion  can  be  measured  using  a  variety  of  techniques  including 
optical  deflection,  optical  interference,  capacitance,  and  tunneling  current.  All  of 
these  techniques  can  be  used  to  produce  an  image  of  the  sample  surface.  As  an 
example.  Figure  2  shows  a  topographical  image  of  nanodots  on  a  polystyrene  surface 
using  an  AFM  with  a  silicon  nitride  tip  [54]. 

Beyond  the  basic  physics  of  a  cantilever  system  and  its  forcing,  there  are  three 
major  methods  to  model  the  forces  that  occur  between  the  AFM  tip  and  the  sample 
surface.  The  first  method  is  to  model  the  tip  surface  interaction  by  means  of  a  po¬ 
tential  function.  A  potential  function  is  a  molecular  model  with  two  terms;  one  term 


Figure  2.  A  sample  topographical  image  of  nanodots  created  on 
a  polystyrene  surface. 


represents  the  adhesive  potential,  and  one  term  represents  the  repulsive  potential  be¬ 
tween  two  molecules.  It  is  a  simple  mathematical  model  that  describes  the  interaction 
force  between  two  molecules  based  on  how  far  apart  they  are.  An  example  of  a  poten¬ 
tial  function  is  the  Lennard- Jones  potential,  which  is  also  called  the  6-12  potential. 
It  can  be  written  as  l/(r)  =  ^  —  ^,  where  r  is  the  distance  between  molecules,  and 
A  and  B  are  constants  based  on  properties  of  the  molecular  interactions  [4], 

The  second  method  is  a  macroscopic  method  that  does  not  explicitly  consider 
atomic  scale  forces.  This  method  is  by  modeling  the  surface  adhesion  between  the 
tip  and  the  surface,  or  between  the  tip  and  the  small  layer  of  condensed  water  on 
the  surface  [78].  The  third  method  is  a  continuum  model  that  uses  the  stress  that 
is  induced  within  the  material  by  the  movement  of  the  AFM  tip  in  contact  with  the 
surface.  This  stress  exerts  a  force  back  onto  the  AFM  tip.  The  effects  of  one  of  these 
methods  to  model  tip  to  surface  interactions,  combined  with  the  properties  of  the 
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surface  material  and  the  basic  physics  of  the  cantilever  system  lead  to  the  ability  to 
measure  surface  material  properties.  Therefore,  we  must  choose  the  most  applicable 
method  to  model  the  interactions  in  order  to  capture  material  properties. 

2.1.1  Cantilever  Beam  Model  With  A  Surface  Potential. 

An  AFM  cantilever  can  be  represented  by  the  shape  of  a  rectangular  beam  with 
a  length  L,  a  width  w  {w  ^  L),  a.  thickness  h  {h  ^  L),  and  a  tip  length  of  Lup. 
We  will  set  up  the  coordinate  system  as  in  Figure  3.  Since  forces  can  occur  in  any 
direction,  the  cantilever  tip  can  deflect  in  any  direction.  We  will  let  the  displacement 
of  the  beam  be  represented  by  u{t,  x,  y,  z),  where  u  =  [u^,  Uy,  u^Y ■  The  terms  Ux,  Uy, 
and  Uz  represent  the  displacement  in  the  x,  y,  and  2;-directions,  respectively. 


Figure  3.  The  model  cantilever  and  the  coordinate  system. 


The  cantilever  tip  movement  is  driven  by  the  forces  that  act  upon  it  and  its 
response  to  these  forces.  In  order  to  study  the  cantilever  tip  movement,  we  must 
first  look  at  the  elastic  properties  of  the  cantilever.  Although  the  cantilever  has 
three  degrees  of  freedom,  for  our  simple  model,  we  will  only  concern  ourselves  with 
the  cantilever  movement  in  the  2;-direction.  Thus,  we  only  need  the  elastic  property 
of  the  stiffness  in  the  ^-direction,  k.  So,  we  will  make  the  following  assumptions. 
First,  we  will  deal  only  with  static  beam  loading.  Next,  we  will  assume  that  there 
is  a  vertical  point  force,  applied  downward  at  the  centroid  of  the  beam  at  the 
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cantilever  tip  end.  This  assumption  means  that  there  is  no  beam  displacement  in 
the  x-direction.  Though  there  is  an  indirect  displacement  in  the  ^/-direction  from  the 
bending  of  the  cantilever  in  response  to  the  vertical  force,  this  displacement  is  of  a 
very  small  magnitude,  and  we  will  assume  that  its  contribution  is  negligible  compared 
to  the  displacement  in  the  2;-direction.  Therefore,  u{t,  x,  y,  z)  can  be  represented  by 
Uz{x,y,z).  The  coordinates  (0,0,0)  represent  the  center  of  the  thickness  and  the 
width  of  the  beam  at  the  clamped  end  of  the  beam. 

Since  the  cantilever  is  only  bending  a  very  small  amount  in  the  vertical  direction 
and  it  behaves  as  a  linear-elastic  material,  the  cantilever  follows  Hooke’s  Law.  This 
states,  Fz  =  kuz,  and  from  Saada  [77],  we  obtain  the  formula 

,  Ewh^ 


where  E  is  the  Young’s  modulus  of  elasticity  for  the  cantilever  material. 

To  ensure  that  our  model  includes  all  of  the  important  physics  of  the  AFM,  we  will 
calculate  the  natural  frequency  of  the  cantilever.  We  will  assume  that  our  cantilever 
is  isotropic  and  has  a  constant  mass  density  with  mass  m.  Drawing  again  from  Saada, 
we  see  that  a  statically-loaded  cantilever  deflected  under  a  vertical  load  at  the  free 
end,  obeys  the  equation 


2F 

M^(0,2/,0)  = -^^(3L-|/)|/2.  (2) 

This  equation  is  for  a  cantilever  under  plane  stress,  and  the  coordinates  (0,  y,  0) 
represent  the  center  of  the  thickness  and  the  width  of  the  beam.  Also,  y  =  0  represents 
the  clamped  end  of  the  beam,  and  y  =  L  represents  the  free  end.  Equation  (2)  assumes 
the  boundary  conditions,  M2(0,0,0)  =  0  and  ^(0,0,0)  =  0.  From  here  on,  we  will 
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let  ^2(0, 1/,  0)  =  Uz{y).  So,  the  cantilever  deflection  at  the  free  end  is 


Uz{L) 


Ewh^' 


(3) 


This  leads  to 


(4) 


If  we  consider  slowly  varying  with  time,  then  the  static  beam  Equation  (4)  captures 
the  dynamic  behavior  as 


Uzit^y) 


(5) 


We  will  next  calculate  the  kinetic  energy,  and  potential  energy,  Ep,  of  the 
cantilever.  First  we  consider  an  infinitesimal  beam  element  dy  a  distance  y  from  the 
clamped  end.  The  infinitesimal  kinetic  energy  of  such  a  beam  element  is  described 
by  the  equation 

dEk  =  ^m{uz{t,y)f^,  (6) 

where  the  dot  is  a  derivative  with  respect  to  time.  Substituting  Equation  (5)  into 
Equation  (6)  and  integrating  over  the  length  of  the  cantilever  produces 


Ek 


1 

2 


m{uz{t,y)Y 


dy 

L 


33  m 

lioT 


(7) 


Since  the  point  force  Fz  acts  only  on  the  free  end  of  the  cantilever,  Ep  is  equal  to  the 
work  needed  to  move  the  free  end  of  the  cantilever  a  distance  Uz{t,  L).  Thus, 


ruz{t,L) 


i'Uz{t,L) 


Ep  = 


Fzduz  = 


kuzduz  = 


kul{t,L) 


(8) 


The  total  energy  of  the  system,  W,  is  equal  to  the  sum  of  the  kinetic  and  potential 
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energies.  If  we  assume  that  the  cantilever  vibrations  occur  without  energy  dissipation 
then  W  is  constant.  If  we  differentiate  W  with  respect  to  time  under  this  assumption, 
and  then  divide  through  by  Uz(i,  L),  we  arrive  at 


33m 

W 


Uz(t,L)  +  kuz(t,L) 


0, 


(9) 


which  is  the  equation  of  movement  for  the  free  end  of  the  cantilever. 

Let  M  =  represent  the  effective  mass  of  the  cantilever.  We  will  suppress  the 
L  and  subscript  ^  symbols  and  let  u{t)  =  Uz{t,L).  Making  the  above  substitutions 
in  Equation  (9)  gives 

Mu{t)  +  ku{t)  =  0.  (10) 


Furthermore,  the  natural  frequency  of  the  cantilever  is  Wq  =  \J 

Next,  we  examine  the  force  of  internal  damping  within  the  cantilever  itself.  The 
internal  damping  is  caused  mainly  by  internal  interatomic  friction.  At  small  velocities 
of  oscillation,  v,  the  force  of  internal  damping  is  modeled  proportional  to  the  velocity 
itself,  i.e.,  Fd  =  —fiov  =  —(3qu,  where  /3o  is  a  positive  constant.  When  the  force  of 
internal  damping  is  included  in  Equation  (10)  we  obtain 


Mil  +  (5qu  +  /cm  =  0. 


(11) 


Dividing  Equation  (11)  by  M  and  make  the  substitution  <5  =  ^,  we  arrive  at 

ii  +  2(5-u  +  cUqM  =  0.  (12) 


The  quality  factor  of  the  system,  or  Q-factor,  is  a  dimensionless  parameter  that 
is  related  to  the  loss  within  the  AFM  cantilever  itself  [8].  It  is  proportional  to  the 
ratio  of  the  stored  energy  of  a  system,  E{t),  versus  the  energy  lost  in  that  system 
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over  a  period  of  time  T.  If  AEt  =  E{t)  —  E{t  +  T),  then  it  is  positive  when  energy 
is  dissipated.  The  Q-factor  characterizes  the  rate  of  energy  transformation  in  the 


cantilever,  and  is  dehned  by 


271  E{t) 

AEj' 


(13) 


When  there  is  very  little  internal  damping  (cjq  ^  5)  the  total  energy  of  the  system 

is 

E(t)  =  Eie-^^\  (14) 


where  Ei 


^McJqZ^  is  the  initial  magnitude  of  the  stored  energy  in  the  system  and 


Z  = 


u(0)  +  (5m(0)  \ 


(15) 


Now,  by  combining  Equations  (13)  and  (14)  and  expanding  the  exponential  function, 
we  get 


271 


277 


.-2(5T 


1  -  (1  -  25T  +  2{6Tf  +  ■  ■  ■ ) 


-  =  -  ^  ^  (16) 

5T  25  25’  ^  ^ 


where  =  ^  =  y/wg  —  5^  and  5T  1.  Therefore,  the  Q-factor  dehnes  the  internal 
damping  of  the  oscillations  in  the  system.  Thus,  Equation  (11)  can  be  rewritten  as 


Mu  + 


Muq 


u  +  ku 


0. 


(17) 


Now,  we  need  to  add  the  oscillations  of  the  cantilever  into  the  model.  There  are 
three  main  methods  to  create  the  cantilever  tapping.  The  hrst  method  is  by  tuning 
the  feedback  control  voltage  on  the  AFM  piezo  crystal  near  the  ringing  conditions. 
The  second  is  by  applying  an  alternating  voltage  to  the  piezo  crystal  in  the  ^-direction. 
The  third  method  is  to  illuminate  the  cantilever  with  a  modulated  laser  at  appro- 
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Figure  4.  The  2-dimensional  cantilever  system  with  a  free  end 
oscillation  of 

priate  frequencies.  For  the  sake  of  simplicity,  all  of  the  preceding  methods  can  be 
mathematically  represented  by  where  B  is  the  amplitude  of  oscillation  and 

u  is  the  frequency  of  oscillation. 

If  we  oscillate  the  free  end  of  the  cantilever,  and  we  translate  the  coordinate  system 
by  shifting  the  beam  in  the  positive  2;-direction  a  distance  uq,  as  shown  in  Figure  4, 
we  get  the  free  end  displacement  equation 

Mil  +  l^ou  +  ku  =  k{uo  +  B  sm{oot)).  (18) 

Next,  we  will  assume  that  the  height  of  the  surface  can  be  represented  by  the 
function  f{y)  as  shown  in  Figure  4.  We  will  also  assume  that  the  surface  has  a 
potential,  V,  that  is  homogeneous  in  the  ^/-direction,  i.e.,  the  potential  only  depends 
on  the  height  above  the  surface  and  not  on  the  location  along  the  surface.  We  are 
adding  a  potential  into  the  model  in  order  to  obtain  a  parameter  to  control  the 
interaction  force  between  the  beam  tip  and  the  surface.  This  force  can  come  from 
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a  number  of  sources  including  molecular  attractions  and  repulsions.  We  note  that 
— Vld  is  the  force  from  the  potential  but,  the  only  pertinent  force  for  our  equation 
is  that  which  is  in  the  ^-direction.  Thus,  adding  the  force  from  this  potential  to 
Equation  (18)  gives 


Mu  +  (5qU  +  ku  =  k{uQ  +  B  sm(ujt))  —  V'{u  —  Uq  —  f{y)), 


(19) 


where  the  prime  means  the  derivative  with  respect  to  the  argument. 

Since  we  are  only  interested  in  the  movement  at  the  tip,  we  will  make  the  substi¬ 
tution,  u{t)  =  y{t)  +  Ltip,  where  Lup  is  the  length  of  the  AFM  tip.  Note,  that  ii  =  y 
and  u  =  y.  Thus,  the  equation  of  tip  motion  is 


My  +  (3oy  +  k{y  Lup)  =  k{uo  B  sm{ut))  -  V'{y{t)  +  Lup  -  Uo  -  f{y)).  (20) 


We  will  now  make  the  assumption  that  our  cantilever  tip  is  always  in  direct  con¬ 
tact  with  the  surface,  which  means  that  uq  =  Lup.  Making  this  substitution  into 
Equation  (20)  and  rearranging  terms  yields 


My  +  Pay  +  ky  =  kB  sm(ut)  -  V'(y{t)  -  /(j/)). 


(21) 


Equation  (21)  is  the  one-dimensional,  simplihed  equation  of  motion  for  the  AFM 
movement  with  a  potential  function. 

The  initial  conditions  are  chosen  so  that  the  tip  has  no  initial  displacement  and  no 
initial  velocity,  i.e.,  |/(0)  =  ?/(0)  =  0.  This  means  the  surface,  /(?/),  is  assumed  to  pass 
through  the  point  {y,z)  =  (0,0).  Now,  divide  Equation  (21)  by  M,  and  substitute 
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/3o  =  ^  and  Wo  =  This  yields 

y  +  ‘^y  +  ^ly  =  sin(wt)  -  ^^\yif)  -  fiy))-  (22) 

Next,  introduce  the  dimensionless  time,  r  =  wot,  and  dehne  the  new  dependent 
variable  x{t)  =  y(t).  Note  that  ij  =  ujqX  and  y  =  oJqX,  where  the  dot  is  differentiation 
with  respect  to  the  argument.  After  making  this  change  of  variables,  and  dividing 
through  by  Wq,  we  obtain 

=  B  sin  (—t]  -  jV'  {x{t)  -  f{y)) .  (23) 

Q  \uJo  J  k 

Further,  the  initial  conditions  become  a:(0)  =  i;(0)  =  0,  and  are  still  quiescent. 

Equation  (23)  with  the  above  initial  conditions  can  be  converted  to  a  nonlinear 
Volterra  integral  equation  of  the  second  kind.  In  Equation  (23)  we  assume  that  Q  >  |, 
which  is  feasible  because  as  we  saw  earlier  in  Equation  (16),  Q  where  6T  -C  1. 

Hence,  our  integral  equation  is 


X(T  = 


/o  yi  -  (1/2Q)2 


F(i(0.0ew(-0siii  (yi-(l/2(J)2(T-0)  d(,  (24) 


where 

F(x(0, 0  =Bsm( -^)  -  V'  (x(0  -  f{y)) .  (25) 

Vwo  J  k 

Equations  (24)  and  (25)  represent  the  cantilever  tip  position  in  terms  of  an  os¬ 
cillating  forcing  function  and  a  potential  function.  We  will  now  proceed  to  look  at 
modeling  the  tip  to  surface  interaction  through  adhesive  forces. 
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Figure  5.  A  diagram  of  the  contact  of  two  spheres. 

2.1.2  Surface  Adhesion. 

Another  method  to  model  the  AFM  tip’s  interaction  with  the  surface  is  through 
adhesion.  Adhesion  is  when  two  surfaces  in  contact  with  each  other,  stick  together 
based  on  the  differing  nature  of  their  materials.  There  have  been  several  different 
surface  adhesion  theories  developed.  The  major  difference  between  these  theories  is 
the  way  in  which  they  idealize  the  contact  interaction  between  two  surfaces. 

Before  contending  with  adhesive  models,  we  will  examine  the  Hertz  model  for  the 
contact  between  two  curved  surfaces.  The  Hertzian  contact  model  [39]  was  developed 
by  Heinrich  Hertz  in  1882.  It  models  the  contact  stresses  and  deformation  that  occur 
when  two  curved  surfaces  are  pressed  together  under  a  constant  loading.  The  model 
assumes  that  the  deformation  of  the  two  surfaces  is  small  compared  to  their  radii  of 
curvature,  and  that  the  materials  of  the  two  surfaces  are  isotropic  and  elastic.  An 
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elastic  isotropic  material  means  that  the  material  has  properties  that  are  identical  in 
any  direction,  and  each  material  can  be  described  by  only  two  elastic  parameters. 

The  Hertzian  contact  model  says  that  the  penetration  depth,  h,  is  equal  to  the 
square  of  the  contact  area  radius,  c^,  divided  by  the  effective  radius,  R  =  .  This 

gives 

h  =  (26) 

where  Ri  and  i?2  are  the  radii  of  each  of  the  spheres.  A  diagram  of  this  is  shown  in 
Figure  5.  Hertz  also  related  the  loading  force,  F,  to  the  contact  area  radius  by 


F  = 


R  ’ 


(27) 


where  k  is  the  effective  Young’s  modulus  of  the  two  materials.  The  effective  Young’s 
modulus,  K,  can  be  calculated  from  the  equation 


1^3  I  ^ 

K  A  El  E2 


(28) 


where  pi  and  p.2  are  the  Poisson’s  ratios  of  each  of  the  spheres,  and  Ei  and  E2  are 
the  Young’s  moduli  of  the  spheres.  The  Hertz  model  is  an  elastic  model  which  does 
not  contain  adhesion.  So,  with  this  model  as  our  basis,  we  will  now  turn  our  focus 
toward  adhesive  models. 

Johnson,  Kendall,  and  Roberts  (JKR)  [48]  used  an  energy  balance  approach  to 
create  the  hrst  model  of  adhesive  forces  between  two  elastic  spheres.  The  force  of 
adhesion,  Fadh,  can  be  written  as 

3 

Fadh  =  ^  -  V^^WaKcl,  (29) 

it 


where  R  =  is  the  effective  radius,  k  is  the  effective  Young’s  modulus,  Wa  is  the 
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work  of  adhesion,  and  Ca  is  the  contact  area  radius.  The  effective  Young’s  modulus,  k, 
is  dehned  by  Equation  (28).  It  is  important  to  note  that  Ca  is  parametrically  dehned 
by  Equation  (29)  and  the  equation 


(30) 


where  h  is  the  penetration  depth. 


The  work  of  adhesion  is  dehned  as 


Wa  =  ll+l2+  712, 


(31) 


where  71  and  72  are  the  surface  energies  of  the  two  adhering  spheres,  and  712  is  the 
interfacial  energy  between  the  two  spheres.  Surface  energy  is  a  material  property 
that  is  dehned  as  the  minimum  energy  per  unit  area  to  form  a  surface  from  the  bulk. 
Higher  surface  energies  indicate  stronger  intermolecular  forces.  The  interfacial  energy 
between  two  solids  is  the  energy  between  the  interface  of  the  two  solids  per  unit  area. 
While  surface  energies  are  easier  to  measure,  the  interfacial  energies  depend  on  what 
two  surfaces  are  adhering  and  are  much  harder  to  obtain.  Therefore,  we  often  estimate 
the  work  of  adhesion.  Berthelot  [80]  used  the  geometric  mean  to  estimate  the  work 
of  adhesion  as 


ir.  =3  -ywuW^, 


(32) 


where  Wn  and  W22  are  the  works  of  cohesion  of  the  two  solids.  The  work  of  cohesion 


is  the  work  per  unit  area  to  break  the  bonds  of  a  solid  with  itself.  So,  Wn  =  271  and 
W22  =  272.  Thus,  Equation  (32)  can  be  rewritten  as 


fha  ~  2,/7i72. 


(33) 
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The  JKR  model  applies  to  tips  with  large  radius  of  curvature  and  cantilevers  of 
small  stiffness.  The  model  works  well  for  highly  adhesive  systems  [34],  The  surface 
forces  involved  within  the  JKR  model  act  only  within  the  contact  area  and  include 
the  influence  of  van  der  Waals’  forces.  The  model  is  based  on  the  assumption  that 
the  cohesive  zone  is  inhnitesimally  small.  The  cohesive  zone  is  the  area  just  outside 
the  region  of  contact  that  is  subjected  to  adhesive  traction. 

Later,  Derjaguin,  Muller,  and  Toporov  (DMT)  [18]  also  solved  the  problem  of  the 
adhesive  forces  between  two  elastic  spheres.  The  force  of  adhesion  for  DMT  can  be 
written  as 

3 

Fadh  =  ^  -  2TTRWa,  (34) 

it 

where  all  the  quantities  are  dehned  the  same  as  in  JKR  except  the  contact  area  radius, 
Ca,  which  is  simply  defined  as 

Ca  =  VhR.  (35) 


The  DMT  model  applies  to  tips  with  small  radius  of  curvature  and  cantilevers  of 
high  stiffness.  The  model  works  well  for  low  adhesion  systems  [34].  The  surface  forces 
involved  within  the  DMT  model  act  only  outside  the  contact  area  and  do  not  include 
the  influence  of  van  der  Waals’  forces. 

Maugis  [62]  then  showed  that  JKR  and  DMT  were  on  opposite  sides  of  the  solution 
spectrum,  and  he  also  came  up  with  an  analytical  solution  to  span  the  transition  from 
JKR  to  DMT.  The  Maugis  solution  can  be  applied  to  any  system  spanning  the  high 
to  low  adhesion  range.  The  amount  of  adhesion  is  determined  by  the  parameter 


^  _  2.06  fRWa^y 


(36) 


where  .^o  is  the  equilibrium  interatomic  distance  between  the  two  spheres.  As  A  — )■  0 
the  system  adhesion  decreases  and  the  stiffness  of  the  materials  increases  so,  the 
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model  behaves  like  DMT.  As  A  — ?•  oo  the  system  adhesion  increases  and  the  materials 
become  more  compliant  so,  the  model  behaves  like  JKR.  Maugis’  model  assumes  that 
the  adhesion  force  acts  within  an  annulus  at  the  contact  area  border. 

Now,  we  will  go  back  to  Equation  (21).  From  here  we  will  insert  the  force  of 
adhesion  from  JKR  in  place  of  the  potential  function.  Since,  we  are  dealing  with  the 
cantilever  tip  radius  in  contact  with  a  surface  which  when  compared  to  the  cantilever 
tip,  can  be  taken  to  be  flat,  i.e.,  R2  =  C)0,  then  R  =  Ri.  Including  the  force  of 
adhesion  in  the  one- dimensional  AFM  model  produces 

3 

My  +  +  ky  =  kB  +  ^  -  ^/67rWaHcl.  (37) 

it 

To  solve  this  equation,  we  will  proceed  as  we  did  in  the  last  section.  First,  we 
are  going  to  choose  the  same  initial  conditions,  i.e.,  y{0)  =  y{0)  =  0.  Next,  we  need 
to  divide  through  by  M  in  Equation  (37)  and  make  the  substitutions  (3o  =  and 

=  \[^-  This  gives 

y  +  ^y  +  ^ly  =  ^Ib  sm{ut)  +  ^  QT^WaKci.  (38) 

Now,  we  introduce  the  dimensionless  time,  r  =  ujQt,  and  dehne  the  new  dependent 
variable  x{t)  =  y{t).  After  making  this  change  of  variables,  and  dividing  through  by 
Wq,  we  arrive  at 

Ht)  +  =  ^sin  f — ^  (39) 

Q  \cJo  J  kR  k 

where  x  =  ^  and  the  initial  conditions  become  x(0)  =  x(0)  =  0. 

Equation  (39)  with  the  above  initial  conditions  converts  to  a  nonlinear  Volterra 
integral  equation  of  the  second  kind.  While  solving  Equation  (39)  we  again  assume 
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that  Q  >  ^,  and  our  integral  equation  emerges  as 


/ 


•r 


where 


=  5sin 


(41) 


We  must  also  remember  that  we  will  need  to  simultaneously  satisfy  the  equation 


h  =  f(v)  - 


(42) 


as  well.  Here,  f{y)  represents  the  height  of  the  surface  at  the  current  cantilever  tip 
position  as  shown  in  Figure  4. 

Equations  (40)- (42)  represent  the  cantilever  tip  position  in  terms  of  an  oscillat¬ 
ing  forcing  function  and  adhesion  defined  by  JKR.  We  will  now  proceed  to  look  at 
modeling  the  tip  to  surface  interaction  through  surface  stresses. 

2.1.3  Surface  Stress  and  Deformation. 

The  final  method  of  producing  a  model  of  the  AFM  is  to  include  the  stress  the 
surface  induces  on  the  AFM  tip  and  the  resulting  surface  deformations.  As  the  AFM 
tip  contacts  the  surface,  it  generates  a  stress,  a,  in  the  surface  based  on  the  magnitude 
and  direction  of  the  force  caused  by  the  tip.  Stress  is  an  applied  force  per  unit  area 
that  causes  a  deformation  per  unit  length,  or  strain  e,  in  a  body.  The  stress  generated 
is  what  causes  the  surface  to  deform.  This  surface  deformation  generates  a  force, 
Fcr{t),  on  the  AFM  tip  as  illustrated  in  Figure  6.  This  leads  to  the  one-dimensional 
AFM  equation 


My  -\-  Poy  +  ky  =  kB  sin  (cut)  -|- 


(43) 
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Figure  6.  Surface  stress  within  the  material,  and  the  resulting 
force  vector  on  the  AFM  tip. 


The  force  F„  includes  the  viscoelastic  surface  forces  as  well  as  the  adhesion  forces. 
This  force  serves  to  couple  the  AFM  model  with  the  material  model,  and  is  highly 
dependent  on  the  geometry  of  the  AFM  tip.  It  can  be  represented  by 


o’zzi'T^O,  z,t)  rdrdO. 

z=0 


(44) 


Here,  azz{T,0,  z,t)  represents  the  surface  stress  normal  to  the  2;-direction  in  cylindri¬ 
cal  coordinates.  The  reason  we  only  include  the  surface  forces  in  the  2;-direction  is 
because  this  is  a  one-dimensional  representation  of  the  AFM  so,  only  the  forces  in 
that  direction  affect  the  model. 

For  the  general  viscoelastic  constitutive  case  [16,  57], 


(45) 


where  Ur,  ug,  and  are  each  functions  of  r,  9,  z,  and  t.  They  represent  the  linear 
displacement  of  the  surface  at  a  given  point,  in  each  of  their  respective  directions. 
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The  displacement  is  written  in  vector  notation  as 


Ur{r,  9,  z,  t) 
u{r,e,z,t)=  ue{r,e,z,t) 
u^{r,9,z,t) 


(46) 


The  Lame  constants,  X{t)  and  /i(t),  represent  the  viscoelastic  properties  of  the  surface. 
Their  time  dependence  means  that  their  past  history  affects  their  current  values. 

2.2  Choosing  the  AFM  Model 

Now  that  we  have  discussed  the  differing  methods  to  represent  the  AFM’s  in¬ 
teraction  with  the  material,  we  can  choose  an  appropriate  model  to  represent  the 
experiment  and  the  material  parameters.  The  AFM  model  is  coupled  to  the  material 
model  by  way  of  the  force  generated  by  the  stresses  within  the  material.  So,  we  will 
generalize  Equation  (43)  a  little  further.  Since,  the  experiment  consists  of  a  static 
preload  and  a  dynamic  forcing,  we  will  introduce  a  static  forcing  displacement,  Bq,  to 
represent  the  displacement  from  the  preload,  and  dehne  B  to  be  the  dynamic  forcing 
amplitude.  In  addition,  we  will  introduce  a  phase  shift,  0,  that  will  be  dehned  later. 
The  phase  shift  occurs  because  the  stress  waves  move  at  differing  rates  in  different 
materials.  So,  the  AFM  model  equation  we  will  use  is 


My  +  (3oy  +  ky  =  k  Bq  +  B  sin  {ut  +  (j))  +  F„{t). 


(47) 


In  summary,  this  chapter  presented  a  mathematical  way  to  represent  the  AFM  as 
a  simple,  spring-mass  system.  We  then  discussed  three  possible  methods  to  model 
the  interaction  of  the  AFM  tip  with  the  sample  surface.  When  this  interaction  was 
modeled  as  either  a  potential  function  or  an  adhesive  force,  an  analytical  solution 
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in  the  form  of  a  Volterra  integral  equation  of  the  second  kind  was  found  for  the 
spring-mass  system.  When  the  interaction  is  modeled  as  a  balance  of  forces  caused 
by  the  stress  in  the  material,  then  we  must  couple  the  AFM  equation  to  a  material 
model  equation.  This  method  allows  the  most  freedom  to  solve  for  particular  model 
parameters  and  to  compare  with  the  experimental  data.  Thus,  the  AFM  model 
coupled  to  a  material  model  will  be  used  to  model  the  experiment. 
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III.  The  Material 


As  a  material  is  acted  upon  by  outside  forces,  structural  changes  and  effects  take 
place  on  and  within  the  material.  These  changes  are  based  upon  the  magnitude  and 
direction  of  the  forces,  as  well  as  the  properties  of  the  material  itself.  This  chapter  will 
present  the  different  models  of  material  behavior  under  load.  Then,  a  general  three- 
dimensional  model  of  a  viscoelastic  material  with  boundary  conditions  is  created  to 
match  with  the  experiment.  The  model  is  non-dimensionalized  to  more  easily  identify 
properties  from,  and  work  with,  the  equations. 

3.1  Material  Behaviors 

Materials  can  demonstrate  three  main  types  of  behavior.  These  behaviors  occur 
under  differing  loading  conditions.  The  material  behaviors  can  be  classified  as  elastic, 
plastic,  and  viscoelastic.  We  call  a  material  an  elastic  material,  a  plastic  material, 
or  a  viscoelastic  material,  if  the  material  exhibits  that  type  of  behavior.  A  material 
may  exhibit  more  than  one  type  of  material  behavior  depending  on  the  loading. 
When  a  load  is  applied  to  a  material  that  behaves  elastically,  the  material  exhibits  an 
immediate  strain.  When  the  load  is  removed,  the  strain  is  immediately  removed  from 
the  material.  We  say  a  material  is  plastic  if  it  exhibits  strain  when  a  stress  is  applied 
but,  when  the  stress  is  removed,  some  of  the  strain  remains  permanently.  Plastic 
deformation  is  defined  to  be  independent  of  time  [24].  The  third  material  behavior 
is  viscoelasticity.  A  viscoelastic  material  exhibits  an  initially  elastic  response  upon 
loading  but,  as  time  goes  on,  there  is  a  slow  and  continuous  increase  of  strain.  When 
the  loading  is  removed,  the  strain  immediately  decreases  but,  not  to  zero.  The  strain 
then  continuously  decreases  from  there  but,  may  never  fully  recover,  leaving  some 
permanent  deformation.  It  is  important  to  note  that  viscoelastic  materials  are  very 
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time  dependent.  The  longer  a  load  is  applied,  the  more  strain  that  will  occur  in  the 
material,  and  the  longer  the  material  will  take  to  recover  after  the  load  is  removed. 

For  this  research,  the  dynamic  nanoindentation  experiment  only  deals  with  poly¬ 
mers  as  the  sample  surface.  Since  polymers  behave  viscoelastically,  we  will  primarily 
concern  ourselves  with  the  modeling  of  these  material  behaviors.  As  the  AFM  tip 
moves  while  in  contact  with  a  viscoelastic  surface,  the  surface  continually  deforms. 
This  deformation  of  the  surface  is  the  primary  cause  of  the  stress  within  the  mate¬ 
rial.  Because  of  the  time  dependent  nature  of  the  material,  and  the  time  dependent 
position  of  the  tip,  the  stress  held  within  the  material  continually  changes. 

In  order  to  ensure  that  we  can  generate  an  analytical  model,  we  will  only  focus 
on  linear  material  models.  Linear  material  models  have  been  shown  to  be  useful 
in  modeling  the  mechanics  of  some  materials  [24].  Many  materials  behave  linearly, 
or  nearly  linearly,  under  small  amounts  of  stress.  A  material  is  considered  to  be 
linear  if,  at  a  given  time,  stress  is  proportional  to  strain,  and  the  principle  of  linear 
superposition  holds. 

Linear  viscoelastic  models  can  be  represented  as  being  made  up  of  two  types  of 
elements,  linear  springs  and/or  linear  viscous  dashpots.  In  these  models,  inertial 
effects  are  neglected.  A  linear  spring  has  the  form 


a{t)  =  Ee{t), 


(48) 


where  E  is  the  Young’s  modulus  constant  of  the  viscoelastic  material.  The  value  of 
E  can  also  be  interpreted  as  a  linear  spring  constant.  A  linear  viscous  dashpot  takes 
the  form 


a(t)  =  ))£((), 


(49) 


where  r]  is  the  constant  coefficient  of  viscosity  of  the  viscoelastic  material,  and  e  is 


the  time  derivative  of  the  strain,  e.  Viscosity  is  the  measure  of  the  resistance  of  a 
fluid  to  deformation. 

Two  simple,  linear,  viscoelastic  material  models,  are  the  Maxwell  model  and  the 
Kelvin- Voigt  model  [24].  The  Maxwell  model  consists  of  two  elements  connected  in 
series.  It  contains  both  a  linear  spring  and  a  linear  viscous  dashpot.  This  model  is 
written  as 

a  a 

The  Maxwell  model  has  no  time  dependent  recovery,  meaning  when  the  stress  is 
removed,  the  linear  spring  element  goes  to  zero  but,  the  linear  viscous  dashpot  re¬ 
mains,  and  is  constant.  This  means  that  there  is  an  instantaneous  recovery  to  a  state 
of  permanent  strain  after  the  load  is  removed. 

The  Kelvin- Voigt  model  consists  of  the  same  two  elements  in  parallel.  It  can  be 
written  as 

.  E  a 
6  -\-  — 6  —  — . 
r]  7] 

This  model  does  not  describe  permanent  strain  after  unloading,  meaning  that  when 
the  load  is  removed,  the  material  eventually  fully  recovers.  This  model  also  does  not 
describe  the  instantaneous  strain  when  a  material  is  loaded  and  the  instantaneous 
partial  relaxation  when  the  load  is  removed. 

Burgers  later  combined  a  Kelvin- Voigt  and  a  Maxwell  model  in  series  to  create 
a  four  element  model  [24].  Zener  also  created  a  similar,  three  element  model,  called 
the  linear  standard  solid  model  [24].  He  did  this  by  combining  a  linear  spring  and  a 
linear  viscous  dashpot  in  parallel,  and  then  he  combined  that  system  in  series  with 
another  linear  spring.  There  are  also  several  other  three  and  four  element  models. 
The  Kelvin- Voigt  and  Maxwell  models,  and  any  other  model  based  on  linear  springs 
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and  linear  viscous  dashpots,  can  be  generalized  into  the  form  [24] 


PqO-  +  pid  +  p2a  H - 


+  Q’lC  +  q2i  +  •  •  •  ; 


(52) 


where  Pi  and  qi  {i  =  0,1,2,.. .)  are  the  appropriate  coefficients  for  the  model  being 
used. 

All  of  these  models  can  be  generalized  even  further  if  we  begin  to  use  fractional 
derivatives.  We  will  start  by  looking  at  the  element 


(53) 


where  p  is  a  proportionality  factor  and  0  <  a  <  1.  When  ck  =  0,  Equation  (53)  refers 


to  a  linear  spring  element,  where  p  is  the  spring  stiffness.  When  a  =  1,  Equation  (53) 
refers  to  a  linear  viscous  dashpot,  where  p  is  the  viscosity.  When  a  is  any  value  in 
between  0  and  1,  we  get  a  transition  element  that  acts  somewhere  in  between  a  linear 
spring  and  a  linear  viscous  dashpot.  Equation  (53)  was  originally  studied  by  Bagley 
and  Torvik  [6].  Koeller  [50]  later  coined  the  term  a  spring-pot.  The  spring-pot  leads 
to  the  general  three  element  fractional  derivative  model,  which  is 


Po(^  +  Pi^(^  =  qo^  +  qi 


(54) 


where  0  <  a  <  1. 

Now  that  we  have  discussed  the  differing  viscoelastic  models,  we  will  be  better 
able  to  choose  a  model  that  best  captures  the  material  properties  we  seek  to  model. 
We  can  now  focus  on  modeling  the  behavior  of  the  material  itself. 
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3.2  Material  Behavior  Model 


Since  the  indentor  of  our  surface  is  the  AFM  tip,  which  is  a  symmetric  geometric 
shape,  we  will  use  cylindrical  coordinates.  The  three-dimensional  equations  from  the 
balance  of  linear  momentum  in  a  viscoelastic  surface  are  [16] 


ft  Q 

[A(f  —  t)  +  2/i(f  —  r)]  —  [VV  •  u]  dr 


d  d'^u 

^u]dT  =  p—,  (55) 


where  p  represents  the  density  of  the  material.  For  a  discussion  on  cylindrical  coor¬ 
dinates,  see  Appendix  A. 

In  addition  to  Equation  (55),  we  need  four  boundary  conditions  for  the  material 
model.  The  hrst  set  of  boundary  conditions  is  at  the  surface  z  =  0,  which  is  where 
the  AFM  interacts  with  the  surface.  At  this  surface,  one  of  the  boundary  conditions 
is  a  mixed  boundary  condition.  That  is,  it  has  a  displacement  representation  from 
0  <  r  <  Ca(t),  and  a  stress  condition  for  r  >  Cait).  Here,  Ca(f)  represents  the  contact 
area  radius  of  the  AFM  with  the  surface.  The  boundary  conditions  a.t  z  =  0  are 


{r,9,0,t)  =  y{t)  -  /(r). 

0  <  r  <  Cait), 

(56a) 

O 

o 

Ca{t)  <  r. 

(56b) 

o 

o 

for  all  r. 

(56c) 

Here,  /(r)  represents  the  geometry  of  the  AFM  tip.  For  the  experiment  we  are 
modeling,  the  AFM  tip  geometry  can  be  approximated  by  a  paraboloid  of  radius  R 
thus,  /(r)  =  1^.  Equation  (56a)  represents  a  coupling  of  the  AFM  to  the  surface, 
and  assumes  that  the  AFM  remains  in  contact  with  the  surface  at  all  times.  It  also 
says  that  the  AFM  tip  acts  as  a  perfectly  rigid  indentor.  Equation  (56b)  represents 
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the  fact  that  outside  the  contact  area  is  a  free  boundary  surface.  If  we  were  to  set  this 
equal  to  some  exponentially  decaying  function,  or  the  negative  gradient  of  a  potential 
function,  it  would  represent  the  force  of  adhesion  for  forces  acting  over  a  short  range, 
and  the  longer  range  forces  (i.e.,  van  der  Waals  forces,  etc.)  acting  at  a  distance. 
Finally,  Equation  (56c)  says  that  the  tangential  surface  tractions  vanish  thus,  /(r) 
must  be  a  smooth  function. 

The  normal  surface  stress,  cT^2(r,  9,  z,  t),  is  defined  in  Equation  (45),  and  the  shear 
stress,  arz{i^,9,  z,t),  is  dehned  as 


arz{r,9,z,t) 


-  t) 


A 

Wr 


dUr  dUz 
dz  dr 


dr. 


(57) 


The  second  set  of  boundary  conditions  occur  at  the  bottom  of  the  polymer  surface, 
where  the  polymer  is  attached  to  an  ideally  stiff  substrate.  The  material  thickness  is 
represented  by  Ml-  Since  the  polymer  material  is  perfectly  bonded  to  the  substrate, 
and  the  substrate  is  perfectly  rigid,  there  is  no  displacement  in  either  the  r  or  2;- 
direction  at  Ml-  This  yields  the  boundary  conditions  at  z  =  Ml  to  be 


r,9,ML,t)  =  0, 

for  all  r. 

(58a) 

r,9,ML,t)  =  0, 

for  all  r. 

(58b) 

Now  that  we  have  setup  the  general  three-dimensional  material  model,  we  can  focus 
on  non-dimensionalizing  both  the  material  model  and  the  AFM  model. 

3.3  Scaling  the  Models 

Before  we  begin  non-dimensionalizing  the  models,  we  must  first  identify  some 
fundamental  properties  of  both  the  AFM  and  the  material.  These  properties  can  be 
seen  in  Table  1.  The  important  AFM  properties  are  the  effective  spring  constant. 
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Table  1.  Typical  parameters  for  the  experiment. 


Parameter 

Value 

R 

30  nm 

k 

0.06  N/m 

UJo 

25  kHz 

Q 

~  100 

p 

1.05  g/cm^ 

Eo 

3  GPa 

1 

V 

3 

k  =  0.06  N/m;  the  resonant  frequency,  tuo  =  25  kHz;  the  quality  factor,  Q  ~  100;  and 
the  tip  curvature  radius,  R  =  30  nm.  The  substrate  material  is  polystyrene,  which 
has  a  density  p  =  1.05  g/cm^,  an  elastic  modulus  of  T'o  =  3  GPa,  and  a  Poisson’s 
ratio  of  =  |.  Now  that  some  of  the  fundamental  properties  are  known,  we  can 
begin  to  scale  the  models. 

In  order  to  scale  the  models,  we  must  look  at  several  natural  length  scales,  and 
several  natural  frequencies,  which  correspond  to  time  scales.  There  are  three  natural 
length  scales  that  arise.  They  are  the  sample  thickness.  Ml]  the  initial  contact  area 
radius,  c°;  and  the  change  in  contact  area  radius,  5a,  caused  by  the  modulation.  Here, 
we  choose  to  write  Cait)  =  c^  +  5ac\{t),  and  note  that  since  the  preload  is  much  greater 
than  the  dynamic  loading,  then  5a  <C  c°.  There  are  also  three  natural  frequencies. 
The  hrst  is  the  undamped  natural  frequency  of  the  AFM,  uq.  This  leads  to  the  time 
scale,  Ti  =  A  =  The  second  is  the  modulation  frequency,  U2.  This  produces 

the  time  scale,  T2  =  Finally,  the  third  time  scale  comes  from  the  elastic  wave 
propagation  time  through  a  distance  i  for  the  elastic  modulus  of  the  surface,  Eq.  This 
leads  to  T3  =  where  ^  will  be  determined  from  one  of  the  length  scales. 

Now,  let  t  =  t*T,  r  =  ^p,  and  2;  =  where  T  and  i  are  to  be  determined 
from  one  of  the  time  and  length  scales,  respectively.  The  time  t*  is  a  dimensionless 
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time,  and  p  and  ^  are  dimensionless  lengths.  It  is  not  necessary  to  scale  6  since 
it  is  already  dimensionless.  However,  we  need  to  scale  our  model  variables  so  they 
are  dimensionless.  To  this  end,  let  y(t)  =  £{qo  +  q(t*)),  Ur{r,6,  z,t)  = 
ue{r,9,  z,t)  =  £v0{p,9,^,t*),  and  Uz{r,9,  z,t)  =  £v^{p,9,^,t*).  The  variables  go  and 
q{t*)  represent  the  dimensionless  static  displacement  and  the  dimensionless  dynamic 
displacement,  respectively.  In  addition,  we  will  scale  out  the  elastic  constant,  Eq,  from 
the  Lame  constants.  Thus,  we  let  p(t)  =  Eopit*)  and  X(t)  =  EoXit*).  Finally,  let 
Cait)  =  c°(l  +  where  (fa  =  %  1  is  an  unknown  dimensionless  parameter, 

and  cl{t)  =  c\{t*T)  = 

Next,  we  will  scale  Equation  (47).  In  addition  to  the  above  scalings,  we  will  let 
/do  =  and  note  that  we  can  write  ^  Also,  we  will  divide  through  by 

to  yield 


m  i  (I)  +  (0 


W*)  +  go] 
+ 


E.{t*T) 

B 

+  —  sin(a;2t*T  +  0) 


(59) 


Now,  applying  these  scalings  to  Equation  (44)  produces 


E.{t*T) 


pdpd9, 

€=o 


where 


(60) 


=  Eq 


[A(r  -  T*)  +  2p{t*  -  r*)] 


+ 


X{t* 


T 


_d_ 

^  dr* 


d 

dr* 

'dVn 


dv^ 


dr* 


_ .  ,  ^  ,  1  dve 

dp  p  p  89 


dr* 


■  (61) 
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Introducing  these  scalings  into  Equation  (55)  and  dividing  through  by  ^  yields 


j‘  [Mt-  +  -T-)]^{VV-v\dT- 

‘  fi(f-T')^lVxVxi7|<iT- 

-CXD 

where  the  cylindrical  coordinate  unit  vectors  are  now  p,  9,  and  The  displacement 
vector  is  now  written  as 


T 


d'^v 


dt 


.2  ’ 


(62) 


v{p,9,i,t*)  = 


ve{p,0,^,t*) 


(63) 


Introducing  the  scalings  into  Equation  (56),  the  boundary  conditions  at  ^  =  0 
become 


v^{p,e,0,t*)  -  qo  + q{t*)  ^  , 

0  <  p  <  y  (1  +  (5 

(64a) 

cr^^(p,6»,0,r)  =  0, 

c 

i 

P  +  6acl{t*))  <  P, 

(64b) 

ap^{p,e,0,t*)  =  0, 

for  all  p. 

(64c) 

where 

/■** 

d 

dvp  dv^ 
dp 

(65) 

(ypi{p,0,i,t*)  =  Eq  1  p{t* 

J  —CO 

—  ) - 

’Ot* 

dr*. 

At  the  bottom  boundary,  ^  Equation  (58)  becomes 

Vp{p,0,  ^,t*)  =  0, 

for  all  p. 

(66a) 

Mp,o,  ^,r)  =  0, 

for  all  p. 

(66b) 

Now,  we  will  choose  i  based  on  our  length  scales.  Since,  we  are  interested  in  what 


35 


occurs  under  the  tip  of  our  AFM,  we  will  choose  £  =  c°.  Next,  we  choose  T  based 
on  our  time  scales.  Before  we  do  so,  we  recall  that  the  experimental  modulation 
frequency,  U2,  occurs  between  1  Hz  and  1  kHz.  Thus,  U2  ujq,  therefore,  Ti  <C 
T2.  This  means  that  the  AFM  will  behave  quasi-statically.  Also,  the  experiment  is 
assumed  to  run  at  a  fixed  frequency  for  a  sufficiently  long  enough  time  so  that  we 
are  at  equilibrium.  So  again,  we  see  that  the  model  will  have  a  quasi-static  solution. 
Looking  at  the  time  scales  in  more  detail,  ~  ~  10“^s,  T2  =  ^  has  a  range  from 
about  10“^s  to  Is,  and  T3  =  c°w  ^  ~  10“^^s.  So,  we’ll  choose  T  =  T2  in  order  to 


capture  the  large  time  scale,  quasi-static  behavior. 

Dehne  A  =  ^  which  ranges  from  25  to  25  x  10^,  and  A  =  ^  = 

which  ranges  from  approximately  10“®  to  10“^^.  Introducing  all  of  the 
above  into  Equations  (59)- (66)  we  get  the  dimensionless  equations  for  the  general, 
three-dimensional  model  of  the  experiment. 

The  dimensionless  equation  for  the  AFM  model  is 


q  +  Q  A^  +  /^i  b  +  qo]  —  Pi  [Bo  +  B  sin(f*  -f-  0)]  xPiFa-, 


(67) 


where  F^(t*)  =  Hq  =  %,  and  H  =  4-  The  parameter  y  =  is  the  ratio 


of  two  forces.  To  see  this,  we  can  write  y  =  •  The  force,  is  proportional 


to  the  elastic  force  the  polymer  exerts  on  the  tip  of  the  AFM  while  the  AFM  is  at 
its  initial  penetration  depth.  The  force,  c^k,  is  the  force  that  the  AFM  spring  would 
exert  when  extended  the  length  of  A  info  fhe  polymer. 

The  dimensionless  equations  for  the  material  behavior  model  are 


(68) 
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and  the  dimensionless  boundary  conditions  are 


and 


n^(p,0,O,r) 

=  %  +  -  f{p), 

0  <  p  <  l+6acl{t*), 

(69a) 

=  0, 

1  +  <  p. 

(69b) 

ap^(p,0,O,r) 

=  0, 

for  all  p. 

(69c) 

Vp{p,9,ML,t*)  =  0, 

for  all  p. 

(70a) 

v^{p,e,ML,t*)  =  0, 

for  all  p. 

(70b) 

where  Ml  =  and  f(p)  =  For  a  parabolic  tip  shape,  fip)  = 

In  summary,  this  chapter  discussed  the  different  types  of  material  behaviors,  and 
the  properties  they  exhibit.  The  main  focus  was  on  the  linear  viscoelastic  behav¬ 
ior  of  materials,  and  how  they  can  be  modeled.  We  created  a  viscoelastic  material 
model  with  boundary  conditions  that  represent  the  experiment.  The  material  model, 
boundary  conditions,  and  the  AFM  model  from  the  last  chapter,  were  then  scaled  to 
capture  the  quasi-static  nature  of  the  experiment  at  the  AFM  tip. 
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IV.  Reduced  Three-Dimensional  Model 


Three-dimensional  material  models  of  physical  systems  are  often  complicated  and 
need  to  be  solved  numerically.  The  goal  of  this  chapter  is  to  make  appropriate  as¬ 
sumptions  so  that  our  model  still  contains  the  pertinent  physics  but,  has  an  analytical 
solution.  The  assumptions  that  we  will  make  for  our  model  are  that  we  can  approxi¬ 
mate  the  material  thickness  by  a  semi-inhnite  half-space,  and  that  we  have  axisym- 
metry.  Before  delving  into  the  math  and  the  reasons  behind  these  assumptions,  we 
will  hrst  present  some  important  background  information. 

4.1  Infinite  and  Semi-Infinite  Half-Space  Models 

Contact  problems  in  the  theory  of  elasticity  generally  assume  that  a  rigid  body 
acts  on  an  elastic  half-space.  More  recently,  it  has  become  increasingly  important 
to  study  viscoelastic  surfaces.  In  a  book  by  Sneddon  [83],  static  theories  based  on 
both  elastic  half-spaces  and  semi-inhnite  half-spaces  are  derived.  N.  N.  Lebedev 
and  la.  S.  Uhiand  [51]  developed  solutions  for  an  axisymmetric  rigid  body  acting 
upon  an  elastic  layer.  They  were  able  to  express  both  the  stresses  and  displace¬ 
ments  in  terms  of  a  single  function  which,  is  the  solution  of  a  Fredholm  integral 
equation.  R.  Y.  S.  Pak  [69]  presented  a  method  based  on  potentials  for  the  dynamic 
response  of  an  elastic  half-space  to  an  arbitrary,  time-harmonic,  hnite,  buried  source. 
An  arbitrary  set  of  transformed  stress-potential  and  displacement-potential  relations 
were  developed  and  can  be  used  in  a  multitude  of  wave  propagation  problems.  He 
specihcally  looked  at  an  embedded  source  of  uniform  distributions.  Later,  Pak  and 
Guzina[70]  derived  fundamental  Green’s  functions  for  the  elastodynamic  half-space 
problem  involving  mixed  boundary  conditions  and  multilayered  media.  They  solved 
three-dimensional  problems  using  a  method  of  displacement  potentials  with  a  set  of 
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generalized  transmission-reflection  matrices  and  internal  source  flelds.  Their  solutions 
encompassed  both  elastic  and  viscoelastic  cases. 

In  the  theory  of  viscoelasticity,  the  contact  problem  is  defined  by  mixed  boundary 
conditions  which  vary  for  the  region  of  contact  and  the  region  of  no  contact.  In 
general,  the  contact  region  varies  with  time.  Lee  and  Radok  [52]  showed  the  solution 
of  the  viscoelastic  counterpart  of  the  Hertz  problem  in  elasticity  can  be  deduced 
from  the  elastic  solution.  The  solution  is  obtained  for  general  linear  viscoelastic 
operators  on  a  half-space.  T.  C.  T.  Ting  [93]  derived  an  explicit  solution  of  an  integral 
equation  which  arises  in  multiple  contact  problems  for  linear  viscoelastic  half-spaces. 
His  solution  is  valid  for  a  contact  area  which  is  a  multiply-connected  region  or  regions. 
Additionally,  the  time-dependent  contact  area  can  be  any  arbitrary  function  of  time. 

4.2  Dual  Integral  Equations 

Dual  integral  equations  often  arise  when  dealing  with  mixed  boundary  value  prob¬ 
lems.  The  general  form  of  dual  integral  equations  is 

Ki{xt)  [1  -|-  w{t)]  ^{t)dt  =  f{x),  0  <  X  <  1, 

K2{xt)^{t)dt  =  g{x),  x  >  1,  (71) 

where  Ki{xt)  and  K2{xt)  are  the  kernel  functions.  These  kernels  can  be  trigonometric 
functions,  Bessel  functions  of  the  first  or  second  kind,  or  modified  Bessel  functions  of 
the  first  or  second  kind.  $(t)  is  the  unknown  function,  and  w{t)  is  an  arbitrary  weight 
function.  Finally,  f{x)  and  g{x)  are  known  functions  that  arise  from  the  boundary 
conditions  of  the  problem  being  solved. 

When  the  problem  being  modeled  has  axisymmetric,  mixed  boundary  conditions, 
such  as  the  experiment  we  are  modeling,  the  dual  integral  equations  that  arise  often 
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have  Bessel  functions  as  their  kernels.  For  this  reason,  we  will  focus  on  dual  integral 
equations  that  have  Bessel  functions  of  the  hrst  kind  as  their  kernel.  The  general 
form  of  dual  integral  equations  involving  Bessel  functions  of  the  hrst  kind  is 

[1  +  w{t)]  ^{t)dt  =  f{x),  0  <  X  <  1, 

Jfj_{xt)^{t)dt  =  g{x),  X  >  1,  (72) 

where  a  and  /3  are  known  constants,  u  and  jj,  are  the  known  orders  of  the  Bessel 
functions,  w{t)  is  an  arbitrary  known  weight  function,  /(x)  and  g{x)  are  known 
functions  valid  over  their  particular  region,  and  <h(t)  is  the  unknown  function  to  be 
determined. 

There  have  been  many  methods  used  to  study  dual  integral  equations  of  this 

class.  Some  methods  are  based  on  particular  solutions  such  as  Titchmarsh  [94], 

where  in  Equation  (72),  w{t)  =  0.  Noble  [65]  used  a  multiplying  factor  method  on 
a  class  of  dual  integral  equations,  where  the  Bessel  functions  are  of  the  same  order 
(z/  =  g).  Nasim  [64]  used  the  properties  of  Mellin  transforms  to  reduce  the  dual 
integral  equations  to  a  single  integral  equation  which  can  be  solved  using  Hankel 
transforms.  Most  methods  use  some  form  of  the  integral  operator  to  reduce  the 
dual  integral  equation  to  a  single  integral  equation,  which  can  be  solved  using  a 
different  method  [61].  Mandal  [60]  used  Sonine’s  integrals  to  reduce  Equation  (72) 
to  a  Fredholm  integral  equation  of  the  second  kind.  Mandal’s  method  is  the  most 
general  and  applies  to  all  the  possibilities  of  Equation  (72).  Since  we  will  use  the 
results  from  Mandal  [60],  a  brief  summary  of  that  paper  is  given  in  Appendix  B. 


40 


4.3  Axisymmetry 


We  will  now  make  assumptions  based  on  the  physics  of  the  experiment  so  that 
we  may  simplify  the  model  for  the  experimental  system.  Recall  the  scaled,  three- 
dimensional  model  of  the  experimental  system  that  we  created  last  chapter.  The  hrst 
assumption  we  will  make  is  that  the  AFM  tip  forces  produce  axisymmetric  loading. 
This  is  a  valid  assumption  since  we  only  have  AFM  tip  movement  in  the  ^-direction, 
and  the  AFM  tip  is  modeled  as  a  paraboloid  of  revolution.  As  a  result,  there  is  no  6 
dependence  for  any  parameter,  ue  is  identically  zero,  and  all  derivatives  with  respect 
to  6  vanish.  Equation  (67)  remains  as 


q  +  +  /3i  b  +  qo]  —  f^i  [Bq  +  B  sin(f*  -1-  0)]  x/5^Fk, 


(73) 


however,  F„  is  now  written  as 


■OO 


pdp. 

C=o 


(74) 


If  we  dehne  the  linear  operator  LpQ  =  ^  +  p  =  (Pfi');  then  Equation  (61)  becomes 


=  Bo 


—  T 


)] 


d  ,  . 

^  vtF  dr 


(75) 


Equation  (68)  remains  as 


(76) 
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where  we  keep  in  mind  that  all  derivatives  with  respect  to  6  vanish,  and  the  displace¬ 
ment  vector  is  now  written  as 


0 


The  boundary  conditions,  Equations  (69)  and  (70),  become 


(77) 


v^{p,0,t*) 

=  Qo  +  Qit*)  -  f{p), 

0  <  p  <  l+6acl{t*), 

(78a) 

a^^{p,0,t*) 

=  0, 

1  +  SaClit*)  <  p. 

(78b) 

=  0, 

for  all  p. 

(78c) 

Vp{p,ML,t*) 

=  0, 

for  all  p. 

(79a) 

V^{p,ML,t*) 

=  0, 

for  all  p. 

(79b) 

respectively.  The  shear  stress  is  now  given  by 


cr 


_d_ 

dr* 


dvp  dv^ 
di  dp 


dr* 


(80) 


Equations  (73)-(80)  represent  the  reduced  axisymmetric  model  of  the  system. 


4.4  The  Hankel  Transform 


We  will  now  make  use  of  the  Hankel  transform  by  dehning 


poo 

Vp{p,^,t*)=  (3Ji{(3p)gp{^,(3,t*)d(3 

Jo 


(81) 
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and 


f5(P.  ?.*■)=  /  PJo(^P)gi(^,l3,(‘)dj3. 


(82) 


The  associated  inverse  Hankel  transforms  of  Equations  (81)  and  (82)  are  defined  to 
be 

poo 

(83) 


and 


9pi^,l3,t*)=  pJi{/3p)vp{p,^,t*)dp 
Jo 


(84) 


The  requirement  for  the  existence  of  Equations  (81)  and  (82)  are  that  the  functions 
\/]3gp{^,  and  fiit*)  be  piecewise  continuous  and  absolutely  integrable 

for  (3  G  (0,  cxo)  [2], 

Applying  Equations  (81)  and  (82),  we  see  that,  Equations  (73)  and  (74)  remain 
unchanged  while.  Equation  (75)  is  transformed  to 


(*t*  noo 


_J  —OQ  ^0 


^JoWp)  [A(t*  -  T*)  +  2p{t*  -  T*)]  r*)dl3dT* 

j-t*  POO  ^ 

+  j  J  /5^</o(/5p)A(f -r*)— ^p({,/3,r*)d^dr* 


,  (85) 


where  '  denotes  and  we  have  used  the  property  LpJi{l3p)  =  [pJi{(3p)]  = 

fJJoifJp)  [12], 

Substituting  Equation  (85)  into  Equation  (74),  interchanging  the  order  of  integra- 
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tion,  and  evaluating  the  p  integral  yields 


F^{t*)=2Tl{l  +  5a^a{t*)) 

ft*  /-oo  _  a 

+  6acl{t*)))  [X{t*  -  T*)  +  2p{t*  -  T*)]  —g'^{0,P,T*)dpdT* 


X 


_J  —  oo  J  0 


+ 


ft*  /-oo  -  (9 

„1  /  +  *  \  \\  \/  .^  -i-X  ^ 


'  —oo  Jo 


fiJim+wn))Ht' -T-)-^g,(o,fi,T-)dpdT- 


.  (86) 


Here,  we  have  taken  advantage  of  the  fact  that  =  0  for  p  >  1  +  6ac\(t*).  We  have 
also  used  the  relationship  pJo{(]p)  =  [pJi(/3p)]  [12], 

Applying  our  Hankel  transforms  to  Equation  (76),  separating  into  two  non-zero 
components,  and  using  the  properties  LpJi{/3p)  =  /3Jo{/3p)  and  ^Jq{(3p)  =  —(5Ji{(5p) 
[12],  produces 


pt*  /-oo  _  fj 

-j  J  PMPp)  -  ^*)  +  -  r*)]  —  (3,  T*)  +  /3,  r*)]  dfddr* 

pt*  poo  o 

+  J  J  /5<7i(/3p)/i(r  -  r*)— [/3p^(c,^,r*) +p"(^,/3,r*)] 

^  Hi  r l}MI}p)^iS,lt,T*)dl3,  (87) 


pt*  /-OO  _  f) 

J  J  PMPp)[\(t'-T-)  +  2li{t'-T‘)]-^,[ldg',(i,P,T-)+g'l{C,P,T-)]dPdT- 

rV  /-CO  o 

-J  J  PM^p)pit* +  Pg'p{^,P,r*)]dpdT* 

=  pU"‘ PMPp)^((,P,T‘)dp.  (88) 

Jo  ot* 

Now,  we  can  take  the  appropriate  inverse  Hankel  transforms  of  Equations  (87)  and 
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(88)  to  obtain 


+ 


[X{t*  -  T*)  +  2fi{t*  -  r*)]  ^  [P'^gpi^,  P,  T*)  +  P,  r*)]  dr* 

I  fl{t* -T*)—[(3g'^{^,f3,T*)  +  g”{^,(3,T*)]dT*  =  (89) 


and 


J  [A(r  -T*)  +  2/i(r  -  r*)]  ^  [fdg'pi^,  (d,  r*)  +  g'l (C,  /3,  r*)]  dr* 

-j  Ht* -T*)-^[(d^g^{^,/3,T*)+/3g'p{^,/3,T*)]dT*  = (90) 

Substituting  Equation  (82)  into  Equation  (78a)  yields 

fdJQWp)g!i{0,/3,T*)d(3  =  qo  +  q{t*)  -  f(p),  0  <  p  <  1  +  dacl{t*),  (91) 

and  substituting  Equation  (85)  into  Equation  (78b),  and  dividing  both  sides  by  Eq, 
yields 

l-t*  1-00  a 

J  J  /dJoildp)  [A(r  -  T*)  +  2p{t*  -  r*)]  T*)dl3dT* 

pt*  _  Q 

+  j  /3^JomX{t*-T*)—gp{0,(d,T*)d/3dT*=0,  l  +  5acl{t*)<p.  (92) 

Applying  the  Hankel  transforms,  Equations  (81)  and  (82),  to  Equation  (80),  and 
using  the  property  ^  Jo(/5p)  =  —(3Ji{(3p)  [12],  admits  the  equation 

pt*  roo  Q 

<^ps(A5.t')  =  Bo  J  J  -  r’)g^  [gpK,/J,r')  -  fejK.A-r')]  ci/So!?-'. 

(93) 

Substituting  Equation  (93)  into  Equation  (78c),  dividing  both  sides  by  Eq,  and  taking 
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the  appropriate  inverse  Hankel  transform  yields 


J  fx{t* -  f3g^{0,(3,T*)]dT*  =0.  (94) 

Since  p  is  a  Lame  parameter,  and  is  a  function  of  the  material  properties,  then  it 
must  not  be  zero  everywhere.  Thus,  Equation  (94)  yields 

g'^{0,f3,T*)-(3g^{0,f3,T*)  =  0.  (95) 

Finally,  applying  our  Hankel  transforms  to  Equation  (79)  we  see  that  it  is  necessary 
for 


gpiML,^,T*)  =  0, 

(96a) 

g^iML,f^,T*)  =  0. 

(96b) 

4.5  The  Fourier  Transform 

In  order  to  solve  Equations  (73)  and  (86)-(96),  we  will  need  to  take  their  Fourier 
transforms.  We  define  the  Fourier  transformation  of  a  function,  f{t*),  as 

/OO 

/(r)e-*-*‘dr,  (97) 

■CO 

and  we  define  the  inverse  Fourier  transform  as 

1  7°°  - 

fin  =  ^J  (98) 

It  is  required  that  fit*)  be  piecewise  smooth  and  absolutely  integrable  on  (— cxo,  cx))  [2]. 
Taking  the  Fourier  transform  of  Equation  (73),  dividing  through  by  and  rear- 


46 


ranging,  yields 


Lo(a;;  Q)g(a;)  =  2'k5{u:)  [Bq  -  go]  +  BS{u:,  0)  +  xF^(a;), 


(99) 


where  d(cj)  is  the  Dirac  delta  function, 


S(u,  =  -iir  [e‘*6(u  -  1)  -  e-*S(u  +  1)]  , 


(100) 


which  is  the  Fourier  transform  of  the  shifted  sine  function,  and 


(101) 


Now,  before  we  take  the  Fourier  transform  of  Equation  (86),  we  will  recall  from 
Section  3.3  that  6a  1.  So,  l  +  SaC^it*)  ~  1,  and  by  a  Taylor  series  expansion  around 
0,  Ji(0(l  +  6acl{t*)))  ^  Ji(0)  +  ^cl{t*)  [Jo{(3)  -  ^2(0)]  +  0{/36af.  Thus,  to  leading 
order,  the  Fourier  transform  of  Equation  (86)  is 

_  _  poo 

Fa{u)  =  2Tr  X{uj)  +  2fl{u)  /  Ji{[6)  +  a{u)/3gp{0,  (6,u)]  d/3  +  0{6a), 

Jo 

(102) 

where  gp{^,(3,u))  and  g^{^,(3,u))  are  the  Fourier  transformations  of  gp{^,(3,t*)  and 
respectively,  and  a{u)  =  In  addition. 


A  (a;)  =  iu  /  X{f)df, 

Jo 


(103) 


and 


fi{uj)  =i0i}  e  jji{t*)dt* . 

Jo 


(104) 


We  dehne  the  Fourier  transformed  viscoelastic  Lame  parameters.  Equations  (103) 


47 


and  (104),  in  this  way  so  that  we  have  the  possibility  of  using  the  elastic- viscoelastic 
correspondence  principle.  As  you  will  see  later,  we  will  not  use  the  correspondence 
principle,  and  in  doing  so,  we  can  attain  an  analytic  solution  to  the  linear  approxi¬ 
mation  of  our  problem. 

Taking  the  Fourier  transforms  of  Equations  (89)  and  (90),  and  collecting  like  terms 
yields 


a{u) 


Q!(a;)  ^ 


(105) 


and 


+  [kl  -  a(a;)/l^]  (3,u;)  +  (3{1  -  a{u;))g'p{^, /3,uj)  =  0,  (106) 


where  a(u))  =  ..  Observe  that  a(u))  =  1  —  2a(u)).  The  dimensionless  shear 

wave  number  is  represented  by  and  the  dimensionless  compressive  wave 

number  is  kl  =  x-  Note  that  kl  =  a(uj)k‘^. 

P  X(iO)-\-2fL(LO)  P  \  /  * 

Next,  the  Fourier  transform  of  Equations  (95)  and  (96),  yield 


^p(0,  /3,  (0,  /3,  w),  (107) 

and 

^p(ML,^,a;)  =  0,  (108a) 

^^(ML,^,a;)  =  0,  (108b) 

respectively.  Again,  taking  advantage  of  the  fact  that  (5a  -C  1  so,  l  +  6acl(t*)  ~  1,  and 
Fourier  transforming  the  linear  approximation  of  the  last  two  boundary  conditions. 
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Equations  (91)  and  (92),  produces  the  dual  integral  equations 


^</o(/^p)^€(0,^,a;)d/3  =  27r(5(a;)  [go  - /(p)]  +  g(^^),  0  <  p  <  1,  (109) 

and 


/3^o(/3p)  [g[(0,/3,a;)  +  Q;(a;)/3gp(0,/3,a;)]  d/3  =  0,  1  <  p,  (110) 

where  we  have  divided  through  by  A(a;)  +  2fi(uj)  in  the  latter  equation. 

Equations  (99)-(110)  represent  the  linear  approximation  to  the  fully  transformed, 
scaled,  and  reduced  model  of  our  experimental  system.  The  equations  show  a  reduced 
three-dimensional  model  of  the  material  behavior  coupled  through  the  boundary  con¬ 
ditions  to  a  one-dimensional  model  of  the  AFM  system.  We  can  now  move  on  to 
solving  this  system  of  equations. 


4.6  Reducing  the  Experimental  Model  to  a  Dual  Integral  Eqnation 


In  this  section,  the  coupled  system  of  second  order  differential  equations  for 
p^(^,  /3,  u)  and  Pp(C,  /3,  oj)  is  reduced  to  a  system  of  first  order  equations.  The  standard 
eigenvalue  and  eigenvector  representation  of  the  solution  is  exhibited  with  expansion 
coefficients  determined  through  the  boundary  conditions.  This  will  lead  to  the  dual 
integral  equation  problem  which  is  evaluated  in  the  next  section. 

First,  Equations  (105)  and  (106)  can  be  transformed  to  the  system 


G' 


0 

0 

1 

0 

1 - 

1 _ 

=  DG  = 

0 

0 

0 

1 

9p 

g'l 

a/3^  -  kp 

0 

0 

1 — 1 

1 

9'i 

Jp  . 

0 

(f-« 

I3{l-a) 

a 

0 

.  9'p  . 

(111) 
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The  solution  to  this  system  is  of  the  form 


4 

G  =  ^  A(/5,  (112) 

i=l 

where  the  Aj’s  are  the  coefficients  determined  from  the  boundary  conditions,  and  A* 
and  Vi  are  the  eigenvalues  and  the  associated  eigenvectors  for  the  matrix  D. 

For  this  system,  dehne 


\  = 

A.  =  (13^  -  kiyi\ 


(113) 


which  represent  the  complex  wave  speeds  associated  with  the  mode  shape  given  by 
Jo  (/dp).  Further,  the  positive  branch  of  the  square  root  is  chosen.  Thus,  the  eigen¬ 
values  are 


A 1  \p , 

A2  =  As, 
A3  =  —As, 
A4  Ap, 


with  associated  eigenvectors 


-1 

-1 

Xp 

Vi  = 

T’ 

\p 

r-/d 

1 

V2  = 

1 - 

Xs  ’ 

r-/d 

-1 

V3  = 

[Ar 

’As’ 

r-i 

1 

Xp 

Va  = 

T’ 

f3 

n  T 


n  T 


(114) 


(115) 
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Inserting  these  valnes,  Eqnations  (114)  and  (115),  into  Eqnation  (112),  and  Qp 
become 


Sp  =  --1  -  Ate’''^]  +  -t  [A^e^A  _  . 

Ap  Ag 


(116) 


To  determine  the  Aj’s,  Equations  (107)-(110)  are  required.  First,  we  examine  the 
boundary  conditions  at  Ml  =  which  are  Equations  (108).  Since  Ml  for  the 
experiment  is  30  nm,  70  nm,  or  220  nm,  and  c°  ~  2  nm  then,  typical  values  of  Ml 
are  in  the  10  -  100  range.  Also  note  that  the  contact  area  radius  is  much  smaller  than 
the  thickness  of  the  sample  thus.  Ml  3>  c°.  Further,  because  the  positive  branch  of 
the  square  root  function  is  used,  3fJ(Ap),  3fJ(As)  >  0,  for  all  (3.  This  means  and 
will  see  exponential  growth  as  ^  becomes  large  while,  and  will  see 

exponential  decay.  Therefore,  to  satisfy  Equations  (108),  both  A2  and  A4  must  be 
exponentially  small.  We  choose  to  set  A2  =  A4  =  0  which,  is  equivalent  to  letting 
Ml  00.  This  is  valid  as  long  as  Ml  ^  c°.  Now,  Equations  (116)  are  reduced  to 


Ary  Ao 


(117) 


Using  Equations  (117)  along  with  Equations  (113)  in  Equation  (107)  produces  the 
relation 


Aj(p,i^)  =  (2y  -  kl)  APP.uj). 


(118) 


The  hnal  coefficient,  A3,  is  found  using  Equations  (109)  and  (110).  To  this  end,  dehne 
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A{P,uj)hy 


-aA{/3,  u)  =  g'^{0,  jS,  u)  +  a/3gp{0,  /3,  uj)  =  -aAg 


(2/3^  -  2/3 

2/3A2Ap  A, 


(119) 


and  w{(3)  by 


[1  +  w{/3)]  A{^,  u)  =  ^g^{0,  /3,  u)  =  -^A^. 


(120) 


In  the  above  definitions  we  have  made  use  of  the  facts  that  a  =  1  —  2a  and  kl  =  a/c^. 

P  * 


Combining  Equations  (119)  and  (120)  yields 


1  +  w{l3)  = 


(121) 


Now,  the  boundary  conditions  given  in  Equations  (109)  and  (110)  are  reduced  to 
solving  the  dual  integral  equations 


f3A{f3,uj)Jo{(3p)d(3  =  0, 


1  <  p, 


(122) 


and 


A{(3,u)Jo{(3p)  [1  +  w{(3)]  d/3 


=  2Tr6{u)  [go  -  fip)]  +q{(^), 


0  <  p  <  1,  (123) 


for  the  unknown  function  A{/3,u).  Having  determined  A{/3,u),  Equations  (118)  and 
(120)  are  used  to  produce 


Ai{(3,uj) 


2/32 


kl 


[l  +  w{/3)\A{/3,uj), 


(124) 
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and 


2  \2 

Aid},  ‘1')  =  — ^  ^ 


(126) 


This  completes  the  solution  to  the  transformed  system,  Equation  (111),  with 
g^{^,(3,u)  and  gp{^,(3,u)  given  by  Equation  (117)  along  with  Equations  (124)  and 
(125).  All  that  remains  is  the  determination  of  A{/3,u)  from  Equations  (122)  and 


(123). 


4.7  Reduction  of  the  Dual  Integral  Equation 

In  this  section,  the  dual  integral  equations.  Equations  (122)  and  (123),  will  be 
transformed  to  a  Fredholm  integral  equation  of  the  second  kind  following  a  method 
presented  by  Mandal  [60]  and  reviewed  in  Appendix  B. 

Comparing  Equations  (122)  and  (123)  to  Equation  (236),  we  see  that  A{(3,u) 
represents  the  unknown  function  <h(/3).  Also,  /(p)  =  /i(p, ca)  =  27iS{u))  [go  —  f{p)]  + 
q{u))  and  g{p)  =  0.  Additionally,  a*  =  0,  (3*  =  —1,  and  p  =  u  =  0.  Thus,  from 
Equations  (238),  A  =  7  =  —  1.  We  choose  s  =  0,  since  that  would  yield  A  +  s  =  —  1  > 
— 1.  So,  from  Equations  (239)  with  s  =  0,  Fi(r)  =  F{r)  and  Gi{r)  =  G{r). 

Now,  from  Equations  (241),  ^  =  £  —  |  and  g  =  m  —  Therefore,  we  choose  1=1 
and  m  =  0  so  that  ^  =  g  =  — |  >  —1.  From  Equations  (240),  G{r)  =  0  and 


(126) 


Recall  that  for  a  parabolic  tip  shape,  /(p) 


Using  this  along  with  /i(p,  cu). 
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Equation  (126)  becomes 


F{r-,uj)=  — — 


vrr  dr 


{27r6{uj)qo  +  q{uj))  /  p{r  -  p  )  ^dp 
Jo 

-  (2vr5(a;)||^  ^  -  p^)"5c?p 


2 

vrr 


27r5(a;)go  +  —  2ti5{uj)  — 


^0^2 


(127) 


Now,  applying  all  of  the  above  information  to  Equation  (237)  yields 


A{p,u)  =  p5 


uo 


rF(r;  uj)J_i{rp)dr 


u^w{u)A{u,  u) 


uJi{u)J_i{p)  —  pJi{p)J_i{u) 


du.  (128) 


If  we  use  the  facts  that  J_i/2(x)  =  ■^/;^cos(x)  and  Ji/2{x)  =  \/:;^sin(x)  then,  we 
can  rewrite  Equation  (128)  as 


A{p^oj)  =  /  \  — F{r-,uj)cos{rp)dr 


TT 


TT 


w{u)A{u,  oj) 


Msin(M)  cos(p)  —  psin(p)  cos(m) 


du.  (129) 


Substituting  Equation  (127)  into  Equation  (129),  and  rearranging  terms  yields 


A{p,uj)  =  - 


27rd{u)  (  -  ^  )  +  q{^) 


sin(p) 


P 


2 

71 


R 


71 


w{u)A{u,  cj) 


p-^  p^ 

Msin(M)  cos(p)  —  psin(p)  cos(m) 


du.  (130) 


In  order  to  simplify  the  notation,  we  will  introduce  the  spherical  Bessel  functions 
of  zero  and  second  order.  These  functions  have  a  relationship  to  the  trigonometric 
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functions  as  follows  [1] 


3o{.x) 

nix) 


sin(a;) 

X 

sin(x)  cos(a:) 

O  ~  O 


sin(a:) 

X 


(131a) 

(131b) 


Therefore,  we  can  rewrite  Equation  (130)  as 


1 

A{p,u:)  =  Co{u:)jo{p)  +  C2{u:)j2{p) - /  w{u)A{u,u:)K{p,u)du,  (132) 

TT  Jo 


where  the  kernel  function  is 


K{p,u)  =  2 


Msin(M)  cos(p)  —  psin(p)  cos(m) 


{u  —  p)(sin(M)  cos(p)  +  sin(p)  cos(m))  +  (m  +  p)(sin(M)  cos(p)  —  sin(p)  cos(m)) 


{u  —  p)  sin(M  +  p)  +  {u  +  p)  sin(M  —  p) 

sin(M  +  p)  sin(M  —  p) 
u  +  p  u  —  p 

=  Jo(m  +  p)  +jo(u-p), 


(133) 


and 


co(a;) 

C2(w) 


2 

71 

2 

71 


27r6{u)  (^qo  -  +  qiu)  , 

2r°l 

2vr<5(a;)^  . 


(134a) 

(134b) 


This  concludes  the  reduction  of  the  dual  integral  equations,  Equations  (122)  and 
(123),  to  a  Fredholm  integral  equation  of  the  second  kind.  Equations  (132)-(134). 
In  the  next  section,  we  will  make  a  linear  approximation  to  the  Fredholm  integral 
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equation  of  the  second  kind,  and  solve  for  the  linear  approximation  of  the  unknown, 

A{p,u). 


4.8  Approximation  of  the  Fredholm  Integral  Equation 

In  order  to  estimate  the  integral  in  Equation  (132),  we  must  hrst  take  a  look 
at  the  weight  function,  tc(/3).  Recall  from  Section  4.5  that  In  addition, 

from  Section  3.3,  we  have  that  /3|  ranges  from  10“^®  to  10“^^.  Since,  uj  and  are 
normalized  quantities,  then  they  are  order  1.  Thus,  ||fcs||  <C  1.  So,  we  will  perform 
an  asymptotic  expansion  for  jS  S>  ||A;s||  on  the  weight  function.  First,  Equation  (113) 
can  be  rewritten  as 


/  ,  2 \  1/2 

/  ,  2  \  1/2 


(135) 


where  we  used  the  fact  that  =  ak^.  Using  a  binomial  series  expansion  with 
(3  ^  ||/cs||  on  Equation  (135)  we  obtain 


a  kl 


a 


k^ 


2/32  8(3^  V/56 


A.  =  /3  1 


(136) 


This  allows  us  to  write  the  asymptotic  expansion  with  f5  /§>  ||fcs||  for  A^Ap  as 


AsAp  — 


(1 


Ikt  f  1  CK^  A  „  f  k^ 


(137) 
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We  may  now  rewrite  Equation  (121)  as 


w{(3)  =  -1  + 


=  -1  + 


=  -1  + 


^  2/32  8  ^4  I  ^6 


4/3^  4M(i_«)_1||(|_q;  +  ^^) +C>(|| 


iM 

i  /34 


fcf 


1 

2  /32  8  /34  “I"  i  /36 


(1 


a) 


ii(i-«+f)  +  o(S 


2(1  -  a)  [  2  ^ 


a  k 


UvJ 

1  _  ^  /  3-2a+a2  \  ^ 

^  /32  4(l-a)  y  (/34 


(138) 


Using  the  geometric  series  expansion  for  ||  <^  1  on  the  last  term  in  Equation  (138) 
admits 


w{f3)  =  -1  + 


2(1 


a) 


akl  f  K 

1 - -  +  O  \  — 

2^^ 


X 


1 

1  j _ 

(5^ 


3  —  2a  + 
4(1  —  a) 


(139) 


Therefore,  we  can  write  the  asymptotic  expansion  of  the  weight  function  for  (3  ^  \\ks 


as 


w{P) 


Wi  + 


M 


3a^  —  4a  +  3 

8(1 -a)2 


where 


Wi 


2a -1 
2(1  -  a)' 


We  will  now  rewrite  Equation  (132)  as 


(140) 


(141) 


1 

A{p,uj)  =  Ao{p,Oi}) - /  wiA{u,uj)K{p,u)du  +  Ii{A),  (142) 

^  Jo 
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where 


Mp,  =  Co(^)Mp)  +  C2(w)i2(p), 


(143) 


and 


Ii(A)  = - /  [w(u)  —  wi]  A(u,u)K(p,u)du. 


TT 


(144) 


We  will  now  put  a  bound  on  Ii{A)  in  order  to  see  if  we  can  approximate  Equa¬ 
tion  (132).  Thus, 


1  p 


[w{u)  —  tci]  A{u,  uj)K{p,  u)du 

\ 

\w{u)  —  tail  \A{u,uj)\  |i^(p,M)|  du 


(145) 


Now,  for  any  point  p,  we  have 


\K{p,u)\  =  |jo(M  +  p)  +  jo(M  -p)|  <  \jo{u  +  p)\  +  |jo(M  -p)|  <1  +  1  =  2.  (146) 


In  addition,  assuming  A{u)  continuous  and  bounded  on  the  interval  u  G  (0,  oo),  then 
let  Ma  =  max„6(o,oo)  |^(w)|-  Thus, 


-^1(44)1  <  — /  \w{u)  —  wi\du 


TT 

=  ‘^Ma 

TT 


Uo 


/3*  /-oo 

\w{u)—wi\du+  /  \w{u)  —  Wi\du 

Jp* 


(147) 


where  (3*  G  (0,  cxd)  is  an  arbitrary  point  to  be  chosen  later. 

Since,  w{u)  and  Wi  continuous  on  ti  G  (0,  /?*),  then  let  =  max„g(o,/3*)  \w{u)  —  Wi 
Therefore, 

/?• 

\w{u)  —  wi\du  <  (148) 
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Also,  for  (5*  3>  ||fcs||  and  from  Equations  (140)  and  (141)  we  obtain 


\w{u)  —  Wi\du  = 


I  /3* 


3q!^  —  4q!  +  3 

8(1 -a)2 


+  o(M 


*  =  (149) 


Hence, 


|/i(A)|  <  -Ma 

71 


=  -MaP* 

TT 


P’M„  +  O  (  P 

+  O  I  ^ 


(150) 


Therefore,  choose  /3*  =  ||A;s||'^,  where  0  <  e  <  1.  This  yields 


Ji(A)|  <  -MaWKW^  [M^  +  O  (||A),||2(1-^))]  . 

71 


(151) 


Thus,  we  can  choose  e  such  that  e  is  as  close  to  1  as  possible  so  that  we  still  satisfy 

(3*  ^  ll^sll-  Hence,  Ii{A)  is  small  and  we  can  approximate  Equation  (132)  to  be 

1 

A{p,u)  =  Aq{p,u) - /  WiA{u,u)K{p,u)du.  (152) 

Jo 

4.9  Solving  the  Approximation  of  the  Fredholm  Integral  Equation 

We  are  now  going  to  solve  Equation  (152)  using  an  integral  equation  Neumann  se¬ 
ries.  For  further  explanation  of  an  integral  equation  Neumann  series,  see  Appendix  C. 
First,  dehne  A*  =  —  and  dehne  the  operator 

roo 

Kf{x)=  /  K{x,u)f{u)du.  (153) 

Jo 

Therefore,  Equation  (152)  is  now  written  as 

A{(3)  =  Ao{f3)  +  \*KA{(3),  (154) 
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where  we  have  suppressed  the  oj  dependence  where  appropriate.  Next,  dehne 


TA  =  Aq  +  \*ka. 

Thus,  Equation  (154)  is  now  written  as 


(155) 


TA  =  A. 


(156) 


Now,  if  01  and  02  are  in  the  subspace  spanned  by  {jo, 42},  then 


11/^01-/^0211  =  11/^(01-02)11  =  11/  //(/3,n)  (0i(m)  -  02(n))  dn||.  (157) 

Jo 


From  Appendix  D,  Equations  (255)  and  (256),  we  obtain 


Thus, 


K{/3,u)  (01  (n)  -  02(m))  du 


7r(0i(0)  -02(0))  • 


(158) 


I|//01  -  /702||  =  IItT  (01  -  02)  II  =  7r||0i  -  02||.  (159) 


Therefore,  we  see  from  Appendix  C  that  for  uniqueness  and  convergence  we  require 
that  |A*|7r  <  1.  Thus,  in  order  to  use  the  integral  equation  Neumann  series  we  require 
that  Itcil  <  1. 

Let’s  assume  for  now  that  |tci|  <  1.  By  the  integral  equation  Neumann  series  we 
obtain 

A  =  Ao  +  X*KAo  +  X*^K^Ao  +  •  •  •  +  X*^K^Ao  +  ■■■  .  (160) 

Using  Equations  (255)  and  (256)  from  Appendix  D,  we  obtain 


KAq  —  ttAq, 


(161) 
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K^Aq  =  tiKAq  =  ti'^Aq, 


(162) 


and 


Thus, 


=  tt^A,. 


(163) 


A  =  Aq-  wiAo  +  {-wif  Ao  H - h  (-wi)'"  Aq^ - =  Ao  ^  (-'«^i)^  •  (164) 


k=0 


Using  the  assumption  that  |wi|  <  1,  the  dehnition  of  an  inhnite  geometric  series,  and 
Equation  (141),  we  can  write  Equation  (164)  as 


1  +  Wi 


=  2(l-«)Ao. 


Here,  we  recall  that  a  =  So,  wi  =  ^7^^  or  =  1 


2(AH-/i) 


(165) 


.  We  will  show  that 


|tai|  <  1  in  Section  4.11. 

With  an  approximation  to  A{/3),  we  can  obtain  all  the  Hj(/3)’s  using  Equations  (124) 
and  (125).  Thus,  using  Equation  (112),  we  can  calculate  G.  Assuming  we  can  take 
the  appropriate  inverses  then,  we  have  a  solution  to  our  reduced,  three-dimensional 
material  model. 


4.10  Finite  Stress  at  the  Contact  Area  Boundary 


In  order  to  relate  the  static  penetration  depth  to  properties  of  the  AFM  and 
surface  system,  we  will  go  back  and  take  a  closer  look  at  the  normal  stress  in  the 
Fourier  transform  domain  at  the  boundary  =  0.  To  do  this,  examine  the  Fourier 
transform  of  Equation  (85)  with  =  0  then. 


%(p,0,a;) 


En 


A(cj)  “h 


PJoiPp)  +  «/3^p(0)]  d(3. 


(166) 
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Next,  recall  from  Equation  (119)  that  ^^(0)  +  a/9^p(0)  =  —aA{(5).  Also,  using  Equa¬ 
tion  (165),  —aA{j3)  ~  — 2a(l  —  a)Ao(/9).  So,  Equation  (166)  becomes 


(J55(p,0,a;)  -2EQfi{uj){l  -  a)  l3Joi/3p)  [cojoi/S)  +  C2j2(^)]  d/3,  (167) 

Jo 

where  we  have  used  the  facts  that  a  =  and  AA/J)  =  Cojo(/9)  +  C2j2(/d)-  Now, 
using  Equations  (257)  and  (258)  from  Appendix  D  yields 


%(p,  0,a;)  -2Eofi{u] 


(1  —  a)  J  Co  -l-  C2  (2  —  3p^) ,  0  <  p  <  1, 

Vi  -  pH  0,  1  <  p. 


(168) 


Substituting  Equations  (134)  into  Equation  (168)  produces 


,  ,  ,  4  ^  ,  (1  —  a) 

(j^5(p,0,u;)  ^ - Eo/i(u;)— == 

TT  L 


2Tr6{u)  (  go  +  ^(1  “  2p^)  )  +  Qi^) 


=  c-|^(PiO)  +  (T^^(p,0,a;), 


(169) 


for  0  <  p  <  1.  Since  the  forcing  function  for  the  experiment  is  made  up  of  a  static 
preload  and  a  dynamic  tapping  piece,  we  would  expect  that  for  this  linear  model,  the 
stress  will  also  be  composed  of  a  static  piece  and  a  dynamic  piece.  Thus, 


- Eofi{u) 

TT 


(1  -  a) 


represents  the  static  component  of  the  stress,  and 


(170) 


»Up,o,oj)  =  -l£:„A(H7)= 

represents  the  dynamic  component  of  the  stress. 

Notice  that  at  the  edge  of  the  contact  region,  p  =  1,  the  stress  becomes  infinite 
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unless  go  is  chosen  appropriately.  In  addition,  the  solution  should  be  continuous  at 
the  contact  boundary.  To  keep  the  static  component  of  stress  hnite  at  the  contact 
boundary,  go  is  chosen  such  that  the  static  stresses  vanish  at  p  =  1.  This  leads  to  the 
choice  go  =  ^-  The  choice  of  hnite  and  continuous  stress  at  the  contact  boundary  is 
the  same  as  the  approximation  made  by  Hertz  [84].  The  equation  for  normal  stress 
at  the  boundary,  ^  =  0,  and  within  the  contact  area,  0  <  p  <  1,  becomes 


-8c° 

ttR 


Eofi{u){l  -  a) 


-  p2  + 


R  q{u) 


p 


(172) 


Although  there  is  no  term  that  we  may  choose  to  balance  the  dynamic  stress 
at  p  =  1,  we  shall  see  that  over  the  total  contact  area,  the  force  generated  by  the 
dynamic  stress  is  hnite.  This  also  means  that  the  dynamic  contact  stress  between  the 
parabolic  tip  and  the  surface,  when  the  static  force  is  much  greater  than  the  dynamic 
force,  is  equivalent  to  that  of  a  rigid  hat  punch  with  the  surface  [90,  84,  91]. 

In  order  to  calculate  the  force  generated  by  the  dynamic  stress,  we  combine  Equa¬ 
tions  (102),  (119),  and  (165)  to  obtain 


/*oo 

F„{oj)  ^  -47rp(a;)(l  -  a)  Ji{(3)  [cojo(/5)  +  C2j2(^)]  d(3,  (173) 

where  we  have  used  the  facts  that  a  =  ^^^1  Ao(/5)  =  Cojo(/9)  +  C2j2{P)-  Now, 

using  Equations  (259)  and  (260)  from  Appendix  D,  and  Equations  (134),  we  can  write 


F„{u)  ^  -  a) 


2TrS(u)  (  ^  I  +  9(1^) 


(174) 


We  have  now  reduced  the  system  of  equations  to  the  solution  of  Equations  (99)- 
(104),  where  the  transformed  force  exerted  by  the  surface  on  the  AFM  tip  is  approx¬ 
imated  as  in  Equation  (174). 
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4.11  Solving  the  Reduced  Model 


Since  we  are  interested  in  solving  for  the  movement  of  the  AFM  tip  in  order  to 
compare  our  mathematical  model  to  the  experimental  data  then,  we  will  solve  the 
reduced  system  of  Equations  (99)-(104)  along  with  Equation  (174).  In  order  to  do 
this,  we  will  need  to  take  the  inverse  Fourier  transform  of  our  system  of  equations.  To 
do  this,  we  need  to  hx  Poisson’s  ratio  to  be  a  constant,  then  we  can  write  A  (a;)  and 
in  terms  of  Poisson’s  ratio,  u,  and  the  Young’s  modulus  of  the  material,  E{u). 
This  also  allows  us  to  calculate  a(u))  =  a  =  where  a  is  no  longer  a  function  of 
frequency. 

For  a  viscoelastic  material,  Poisson’s  ratio  is  usually  a  function  of  frequency  be¬ 
cause  the  transverse  strain  may  be  out  of  phase  with  the  longitudinal  strain.  For  the 
solution  of  the  current  system  of  equations,  we  assume  that  this  difference  in  strains 
is  minimal  so  that  Poisson’s  ratio  is  a  constant  with  respect  to  frequency.  It  will  also 
be  important  to  note  that  —l<u<  .5.  So,  from  Findley  et  ah  [24]  we  have 


A(a;) 


—uE{uj) 

(l  +  i^)(2z/-l) 


(176) 


and 


/i(u) 


2(1  +  1^)' 


(176) 


This  admits  to 


1  -  2z/ 

^  2(1-0' 


(177) 


At  this  point,  we  shall  take  a  look  at  wi  =  2{\+fi) '  f^^^all  from  Section  4.9  that 
we  require  |tci|  <  1.  Using  Equations  (175)  and  (176)  we  obtain 


Wi  =  —v. 


(178) 
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Thus,  |wi|  <  1. 


Substituting  Equations  (175)-(177)  into  Equation  (174)  yields 


-2E{u) 

1  —  z/2 


2ti5{uj) 


+  ?(u) 


(179) 


The  static  piece  from  the  above  equation  is  precisely  the  viscoelastic  correspondence 
of  the  Hertz  solution  in  terms  of  force  [84], 

Now,  substituting  Equation  (179)  into  Equation  (99),  using  the  fact  that  go  = 
and  collecting  terms  gives 


[Lo{u)  +  Po(<^)]  q{(^)  =  2tt6{uj) 


/  2 


BS{u,(t)),  (180) 


where 


Po(a;) 


2E{uj) 

(T^^- 


(181) 


Recall  from  Section  3.3  that  y  = 

So,  taking  the  inverse  Fourier  transform  of  Equation  (180),  and  equating  the  static 
and  dynamic  terms  gives 


1  +  3fo(0) 


(182) 


and 


qin  = 

Equations 

system. 


B 


2i(Lo(l)  +  -Po(l)) 


B 


2*(To(— 1)  +  Po{—l)) 


(181)-(183)  represent  the  solution  to  the  scaled,  reduced. 


three-dimensional 
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4.12  Analysis  of  the  Reduced  Three-Dimensional  Model 

Let’s  start  off  by  taking  a  closer  look  at  L^iu)  from  Equation  (101).  First,  we 
recall  that  /3i  is  the  ratio  of  the  resonant  frequency  of  the  AFM,  cjq,  to  the  modulation 
frequency,  a;2,  so  that  /3i  =  ^  is  in  the  range  of  25  to  25  x  10^  which,  is  large  compared 
to  unity.  Also  note  that  Q  ~  100.  So,  Equation  (101)  can  be  approximated  by 
Lq{!jj)  1  for  o;  =  0,  ±1. 

In  addition,  dehne  which  is  the  elastic,  plane  strain  Young’s  modulus. 

Thus,  Equation  (181)  can  be  written  as 


Pq{oj)  =  2x*E{oj), 


(184) 


c'iE* 


where  y*  = 


is  a  dimensionless  ratio.  Therefore,  we  can  write  Equation  (182)  as 


B,  =  ^  + 


R  3Rk 


.£(0). 


(185) 


Here  we  note  that  £"(0)  represents  the  static,  normalized  viscoelastic  Young’s  modulus 
so,  E(0)  =  1.  Thus,  Equation  (185)  becomes 


Bo  = 


^2 

R  3Rk 


(186) 


So,  the  difference  between  Bo  and  go 


c° 

Z2l 

R 


yields 


B 


R 


E* 

-^0 

3Rk 


(187) 


Here,  go  =  ^  represents  the  dimensionless  penetration  depth  of  the  AFM  tip  into  the 
surface,  and  Bq  represents  the  dimensionless  distance  caused  by  the  static  loading, 
i.e.,  the  preload.  Thus,  the  left  hand  side  of  Equation  (187)  represents  the  difference 
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of  the  static  loading  and  the  tip  penetration  into  the  material.  We  know  that  this 
difference  is  caused  by  various  physical  properties  of  the  AFM  and  the  material. 
Examples  of  these  are  the  AFM’s  tip  curvature  radius  and  its  effective  stiffness,  as 
well  as  the  contact  area  radius  and  the  Young’s  modulus  of  the  material. 

The  right  hand  side  of  Equation  (187),  is  a  ratio  of  forces.  The  numerator  rep¬ 
resents  the  force  the  material  exerts  on  the  tip  over  the  whole  contact  area.  The 
denominator  represents  the  force  of  the  effective  AFM  spring  constant  the  length  of 
the  tip  curvature  radius.  In  other  words,  we  have  the  ratio  of  the  force  the  surface 
exerts  on  the  tip  to  the  force  the  tip  exerts  on  the  surface.  Note,  that  all  the  terms 
in  this  ratio  are  properties  of  the  AFM  and  the  material,  as  we  expected. 

Moving  on  to  the  dynamic  portion  of  the  system,  we  see  that  Equation  (183) 
becomes 

_ _  I  ghi*+0)  _  I  ^  ^ 

E{l)j  \l  +  2x*E(-l) 

If  E{1)  and  E{—1)  are  complex  conjugates,  which  is  the  case  for  the  viscoelastic 
Kelvin- Voigt  and  Maxwell  models,  as  well  as  the  generalized  viscoelastic  polynomial 
model,  then  we  can  write  Equation  (188)  as 


l  +  2x* 


(188) 


qin  = 


B 


\l  +  2x*E{l)\ 


sm{t*  +  (j)  +  0i), 


(189) 


where 


A 


h  =  arctan 


1 


i+2x*£;(i)^ 


3f? 


1 


(190) 


.i+2x*£;(i)^ 

Here,  we  see  that  the  amplitude  of  the  dimensionless,  dynamic  penetration  depth, 
q(t*),  is  directly  proportional  to  the  dynamic  forcing  amplitude,  B.  We  also  see 
that  the  amplitude  depends  on  material  properties  and  AFM  properties  through  y* 
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and  E{1).  Both  of  these  observations  make  sense  physically.  We  would  expect  the 
dynamic  penetration  depth  to  be  smaller  in  magnitude  than  the  dynamic  forcing 
amplitude  because  the  surface  will  resist  the  tip,  causing  the  tip  to  not  penetrate 
to  the  full  depth  of  the  forcing  amplitude.  We  would  also  expect  that  the  material 
properties  and  the  AFM  properties  would  help  to  determine  the  penetration  depth. 

We  also  see  that  the  dynamic  tip  movement  is  sinusoidal  as  is  the  forcing  but,  the 
material  properties  and  the  AFM  properties  add  a  phase  shift.  This  is  to  be  expected 
as  the  system  is  not  an  instantaneous  elastic  system  but,  rather  it’s  a  viscoelastic 
system.  We  also  observe  that  (pi  is  a  function  of  our  forcing  frequency,  U2,  by  means 
of  the  viscoelastic  modulus,  E{1). 

This  chapter  presented  some  background  material  on  inhnite  and  semi-infinite, 
three-dimensional  material  models  and  dual  integral  equations.  We  then  made  ap¬ 
propriate  assumptions  in  order  to  reduce  the  system  of  equations  that  represents  the 
experiment.  These  assumptions  eventually  allowed  us  to  solve  the  equations  in  the 
transform  domain.  We  then  inverted  the  transformed  equations  to  obtain  the  solution 
to  the  system  in  an  analytical  form.  Finally,  we  analyzed  the  system.  In  the  next 
chapter,  we  will  present  a  new  model  that  will  allow  us  to  bring  together  the  system 
model  and  the  actual  experimental  setup.  We  will  then  compare  the  combined  model 
to  experimental  data. 


V.  Application  of  Models  to  AFM  Measurements 


In  order  to  compare  the  reduced  three-dimensional  model  to  the  AFM  experimen¬ 
tal  data,  we  will  introduce  a  model  called  the  “error  model”.  This  model  will  allow 
us  to  solve  for  the  penetration  depth  of  the  AFM  tip  into  a  material  using  AFM 
experimental  data.  This  will  allow  us  to  directly  equate  the  analytical  system  model 
solution  for  the  penetration  depth  to  the  error  model  solution  for  the  penetration 
depth.  This  will  yield  an  analytical  solution  for  the  viscoelastic  material  parameters 
in  terms  of  AFM  experimental  data.  We  note  that  the  analytical  solution  may  have 
wide  applicability  to  other  problems  but,  we  will  make  application  to  a  specific  prob¬ 
lem.  In  particular,  we  will  present  results  from  data  collected  for  a  polystyrene  thin 
him  on  a  silicon  substrate. 

5.1  Error  Model 

Recall  that  the  experiment  we  have  been  modeling  is  described  in  section  1.2.  Let 
the  initial  state  of  the  system  be  with  the  AFM  cantilever  tip  ‘just  resting’  on  the 
sample  surface,  i.e.,  with  no  applied  loading.  Let  the  surface  to  air  boundary  be  at 
z  =  0,  with  the  positive  2;-direction  in  the  inward  normal  direction  to  the  surface.  A 
laser  beam  is  rehected  off  the  tip  end  of  the  cantilever  and  onto  the  AFM  detector. 
In  the  initial  state,  the  AFM  detector  records  a  small  voltage  due  to  surface  adhesion 
and  other  small  forces  between  the  surface  and  the  AFM  tip.  The  AFM  detector 
actually  measures  the  position  of  the  top  of  the  cantilever.  However,  assuming  the 
length  of  the  tip  of  the  cantilever  remains  constant  throughout  the  experiment,  the 
voltage  recorded  at  the  AFM  detector  is  the  position  of  the  tip  relative  to  the  original 
surface,  denoted  z(t-,  U2),  where  t  is  time  and  u)2  is  a  tapping  frequency. 

In  the  experiment,  an  initial  DC  voltage,  Qdc,  is  applied  to  the  piezo.  Assuming 
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the  silicon  base,  which  rests  atop  the  piezo,  is  perfectly  rigid  and  perfectly  bonded 
to  the  sample,  the  piezo  movement  is  transferred  into  the  sample.  Further,  the  top 
surface  of  the  sample  moves  the  exact  same  distance.  This  distance  is  determined  via 
through  the  calibration  constant  Qp,  which  converts  the  voltage  of  the  piezo  into 
a  displacement.  For  the  piezo  used  in  this  experiment,  gp  =  g||g  volts  per  nanometer. 
The  AFM  bends  against  and  slightly  penetrates  the  surface  but,  because  of  the  vis¬ 
coelasticity  of  the  surface,  there  is  some  creep  after  the  initial  voltage  displacement. 
After  some  time,  the  AFM  tip  and  the  surface  reach  “near”  equilibrium,  i.e.,  steady 
state  with  a  small  drift.  When  this  equilibrium  position  is  measured  at  the  AFM 
detector,  we  call  this  the  set  point  voltage,  or  preload  voltage,  L. 

A  manual  gain  is  set  at  the  AFM  detector  called  the  modulation  voltage  bias,  7. 
This  tells  the  closed-loop,  computer  controlled  feedback  system  within  what  range  of 
the  preload  voltage  to  achieve,  when  measured  at  the  AFM  detector.  After  equilib¬ 
rium  is  reached,  an  initial  sinusoidal  voltage  is  applied  to  the  piezo  at  a  frequency  of 
UJ2.  The  closed-loop,  computer  controlled  feedback  system  determines  an  amplitude 
of  forcing  and  a  phase  shift  that  modifies  the  initial  modulation  in  order  to  maintain 
the  set  point  voltage  to  within  the  modulation  voltage  bias.  The  modulation  volt¬ 
age  bias  is  established  as  the  reference  phase.  The  first  lock-in-ampliher  records  the 
AC  error  voltage  amplitude,  P{u)2),  and  phase,  9p{u)2),  which  are  the  differences  in 
voltage  amplitudes  and  phases  between  the  modulation  voltage  bias  and  the  voltage 
measured  at  the  AFM  detector.  This  represents  the  error  signal  between  the  can¬ 
tilever  tip  location,  as  driven  by  the  computer  controller,  and  the  actual  tip  position. 
This  error  signal  is  attributable  to  properties  of  the  surface  and  the  piezo  system. 

In  addition  to  the  calibration  constant  gp,  which  converts  piezo  voltages  into  dis¬ 
placements,  there  is  another  calibration  constant  ga,  which  converts  the  voltages  of  the 
error  signal  and  modulation  voltage  to  cantilever  displacements.  For  the  experiment 
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being  examined  ga  =  ^  volts  per  nanometer. 

The  DC  part  of  the  error  displacement,  is  the  difference  between  the  set 

point  displacement,  and  the  static  AFM,  z{t;0),  i.e., 

^  =  --z{t-,0).  (191) 

9a  9a 

At  equilibrium,  =  0,  so  it  follows  that  z{t;  0)  = 

The  entire  experimental  system  is  nonlinear,  which  means  the  total  response  of 
the  system  has  multiple  natural  modes  of  oscillation,  i.e.,  the  tapping  frequency  is 
the  sole  forcing  function  whereas  the  resultant  system  response  occurs  at  multiple 
frequencies,  including  the  tapping  frequency.  For  the  experiment,  the  response  of  the 
error  and  the  computer  controlled  feedback  is  measured  by  means  of  a  LIA  only  at 
the  forced  oscillation  frequency,  U2-  Thus,  the  error  of  the  system  can  be  represented 
by  ^  ^ 

Pp£  0p{u2))  =  —  +  —  sin  {u2t)  -  z{t-,  0J2).  (192) 

9a  9a  9a  9a 

The  left  hand  side  of  Equation  (192)  represents  the  error,  while  the  first  two  terms 
on  the  right  hand  side  represent  the  piezo  movement.  Finally,  z{t]U2)  is  the  AFM 
movement  relative  to  the  original  surface  position. 

Solving  Equation  (192)  for  the  position  of  the  cantilever  tip  produces 


where 


and 


z{t) 


- ^ - sin  {uj2t  +  9c) , 

9a  9a 


C  = 


^2  p2  _  2^p  (.Qg 


1/2 


6  c  =  arctan 


-Psin(6>p)  A 
7  -  Pcos(6>p)  J 


(193) 


(194) 


(195) 
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Note,  we  have  suppressed  the  0J2  dependence  on  the  appropriate  variables,  and  taken 
advantage  of  the  fact  that  Pdc  =  0.  The  term  C  represents  the  voltage  amplitude  of 
the  AFM  tip  position,  and  9c  represents  the  phase  shift  of  the  AFM  tip  in  reference 
to  the  driving  signal.  Both  terms  can  be  calculated  from  experimental  data. 

In  order  to  relate  the  reduced  three-dimensional  model  from  Chapter  IV,  to  the 
error  model,  it  is  necessary  to  relate  the  AFM  tip  movement  to  the  penetration  depth. 
To  this  end,  let  z{t)  =  Zo  +  Zi(t)  be  the  AFM  tip  displacement,  and  y(t)  =  yo  +  yi(t) 
be  the  penetration  depth,  then  zo  and  yo  represent  the  static  tip  location  and  the 
penetration  depth  of  the  cantilever  tip  based  solely  on  the  preload  voltage,  respec¬ 
tively.  Further,  Zi{t)  and  yi{t)  represent  the  dynamic  portion  of  the  cantilever  tip 
movement  and  the  dynamic  penetration  depth  based  on  the  modulation  signal,  re¬ 
spectively.  When  the  system  reaches  dynamic  equilibrium,  the  surface  of  the  material 
in  contact  with  the  AFM  has  moved  in  the  negative  ^-direction,  a  distance  So, 
the  static  penetration  depth  of  the  AFM  tip  is 

Qdc  Qdc  L 

yo  = - Zo= - . 

9p  9p  9a 

Here,  yo  represents  the  positive,  static,  penetration  depth  of  the  AFM  tip  into  the 
sample  surface. 

After  static  equilibrium  is  reached,  the  initial  modulation  and  the  computer  feed¬ 
back  voltages  are  applied  to  the  piezo.  The  second  lock-in-ampli£er  measures  this 
combined  piezo  forcing.  Specifically,  the  LIA  measures  the  piezo  amplitude  voltage, 
Q{u2),  and  piezo  phase  shift,  9q{u2).  Since  the  piezo  is  a  crystal,  there  is  a  time 
delay  associated  with  contraction  or  expansion  due  to  an  applied  voltage.  This  delay 
in  contraction  or  expansion  is  called  the  piezo  delay  time,  Tp{uj2),  and  it  is  a  func¬ 
tion  of  the  forcing  frequency.  Because  of  the  assumed  perfectly  rigid  silicon  base, 
and  perfect  bonding  to  the  sample,  all  the  piezo  movement  is  transferred  into  surface 
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movement.  This  surface  movement  doesn’t  occur  until  after  the  piezo  delay  time. 
Thus,  the  dynamic  AFM  tip  position  is  not  affected  until  after  the  piezo  delay  time. 
This  leads  to  the  equation  for  the  dynamic  penetration  depth  as 


sin  {u2t  +  9c)  ■  (197) 


Here,  yi(t)  is  the  oscillation  of  the  tip  penetration  depth  around  the  static  tip  pene¬ 
tration  depth,  hq. 

Recalling  from  Section  3.3  that  y(t)  =  c°(go  +  Q(t*)),  Equations  (196)  and  (197) 
can  be  written  as  dimensionless  equations.  This  leads  to  the  cantilever  movement 
based  on  the  error  model  as 


(198) 


<?o  —  Qdc  —  E, 


from  Equation  (196),  and 


q(t*)  =  Q  sin  [t*  -  T*  +  9q)  -  C  sin  (t*  9c) , 


(199) 


T*  =  072 Tp.  The  terms  Qdc  and  L  are  the  dimensionless  static  distances  generated 
by  the  piezo  and  the  preload,  respectively.  The  terms  Q  and  C  are  the  dimensionless 
dynamic  displacement  amplitudes  generated  by  the  piezo  and  the  AFM  movement, 
respectively.  Finally,  t*  is  dimensionless  time,  and  t*  is  the  dimensionless  piezo  delay 
time. 
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5.2  Equating  the  Models 


Now,  the  dimensionless  penetration  depth  from  the  reduced  three-dimensional 
model,  and  the  dimensionless  penetration  depth  from  the  error  model,  can  be  equated. 
Before  doing  so,  we  must  hrst  dehne  i?o,  B,  and  0,  as  in  Equation  (67),  in  terms  of 
the  error  model  parameters.  Recalling  Equation  (67),  we  see  that  Bq  represents  the 
dimensionless  static  loading  that  occurs  at  the  AFM,  B  represents  the  dimensionless 
dynamic  amplitude  that  occurs  at  the  AFM,  and  0  represents  the  phase  shift  at  the 
AFM,  as  compared  to  the  reference  phase.  Since,  Qdc  represents  the  dimensionless 
displacement  due  to  static  loading  at  the  piezo,  and  the  substrate  is  assumed  perfectly 
bonded  to  a  perfectly  rigid  silicon  base,  the  movement  at  the  top  {z  =  0)  of  the 
material  surface  is  Qdc-  So,  the  dimensionless  static  forcing  is  Bq  =  Qdc-  By  the 
same  argument,  the  dimensionless  dynamic  amplitude  of  the  system  is  R  =  Q.  When 
a  voltage  is  applied  to  the  piezo  crystal,  movement  at  the  top  surface  of  the  sample 
is  delayed  by  the  piezo  delay  time,  Tp.  In  addition,  the  voltage  applied  to  the  piezo  is 
phase  shifted  from  the  reference  signal  by  9q.  Thus,  the  phase  shift  is  (j)  =  9q  —  t*. 

The  static  portion  of  the  reduced  three-dimensional  model.  Equation  (186),  is  now 
written  as 


Qdc 


cl  4cf-E; 
R  3Rk 


(200) 


Substituting  Equation  (200)  into  Equation  (198),  and  using  the  fact  that  Qq  =  we 
can  solve  for  the  dimensionless  preload  as 


L 


4cfRo*  4  , 


(201) 


Equation  (201)  says  that  the  dimensionless  preload,  L,  is  equal  to  a  constant  times 
the  ratio  of  the  force  the  surface  exerts  on  the  tip,  c°^Rq,  to  the  force  the  tip  exerts 
on  the  surface,  Rk.  Using  Equations  (198)  and  (201),  we  can  solve  for  the  unknown. 
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X*,  producing 


3 

L 

3 

L 

4 

_Qdc  ~  L 

4 

(202) 


Equation  (202)  allows  for  the  calculation  of  x*  based  on  experimental  data.  This 
constant,  x*;  is  independent  of  frequency  and  is  strictly  a  positive,  real  number. 
Additionally,  using  =  ^,  we  can  solve  Equation  (198)  for  c°  to  obtain 


c 


0 

a 


(203) 


Using  Equations  (202)  and  (203),  and  solving  for  Eq  based  on  the  experimental  data 
results  in 


E*o  = 


3kL 


'9a 

1/2 

—Qdc  ~  L 

-R. 

.9p 

1  -3/2 


(204) 


The  dynamic  portion  of  the  reduced  three-dimensional  model.  Equation  (189),  is 
now  written  as 

^  A - +  ^Q-  (205) 

|l  +  2x*E(l)| 

Substituting  Equation  (205)  into  Equation  (199)  yields 


|1  +  2x*E(l)|  +  <^i)  =  <5  sin  {t*  +  0q-t*)  -  Usin  (U  +  dc)  ■  (206) 

Here,  the  left  hand  side  of  Equation  (206)  represents  the  dynamic  portion  of  the 
penetration  depth  from  the  reduced  three-dimensional  model,  and  the  right  hand 
side  represents  the  dynamic  portion  of  the  penetration  depth  from  the  error  model. 
For  convenience,  let  x  =  t*  +  9q  —  r*,  and  combine  the  sine  terms  on  the  right  hand 
side  of  Equation  (206),  to  obtain 


l  +  2x*E(l)| 


sin(x  -|-  0i)  =  D  sin  (x  -|-  ^2) , 


(207) 
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where 


D=[Q‘^  +  C^-  2CQ  cos  {Oc  -eQ  +  r;)] , 


(208) 


and 


02  =  arctan 


/  -Csin^ec -Oq  +  t*)  \ 
VQ  -  Ccos(6'c  -Oq  +  t*)}  ' 


(209) 


For  equality  to  hold  in  Equation  (207),  the  amplitudes  and  the  phases  must  match. 


This  leads  to 


Q 

l  +  2x*E(l)| 


(210) 


and 


01  =  arctan 


1 


l+2x*E(l), 


1 


l+2x*E(l), 


=  arctan 


/  -Csinjec -Oq  +  t*)  \ 
VQ  -  Ccos(6»c  -0q  +  t*)) 


02.  (211) 


From  Equation  (211)  we  get 


(i+2x*£:(i))  _  —Csm.{6c  —  0Q  +  T*) 

Q  -  C cos{ec  -  9q  +  T^y 

\l+2x*E[l)) 


(212) 


Equations  (202),  (210),  and  (212)  represent  the  analytical  solution  to  the  experimental 
system.  The  remaining  unknowns  to  be  solved  for  are  the  model  parameters  from 
whichever  viscoelastic  model  we  choose  to  represent  the  sample  surface. 


5.3  Analysis  of  the  Experimental  Model 

Now  that  we  have  the  solution  to  the  experimental  system.  Equations  (202),  (210), 
and  (212),  we  can  proceed  to  analyze  the  model  and  solve  for  the  unknown  viscoelastic 
material  parameters.  First,  observe  from  Equation  (202),  we  have  the  solution  for  y* 
in  terms  of  the  experimental  data,  L  and  Qdc-  Next,  we  analyze  Equations  (210)  and 
(212),  by  letting  the  term  1  +  2y*E(l)  =  a  +  60  a  complex  number.  Note,  that  both 
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a  and  b  can  be  freqnency  dependent,  and  a  =  1  +  2x*3?(^(l))  and  b  =  2x*‘^{E{l)). 
Thus,  we  can  write  Equation  (212)  as 


-C  sin(6'3: 


a  Q  —  C  cos{9^)  ’ 


(213) 


where  we  have  let  9^  =  9c  —  Oq  +  t*  Solving  this  equation  for  b  in  terms  of  a  yields 


b  =  a 


C  sin(9x 


Q  —  C  cos{9j;) ' 


(214) 


Turning  our  attention  to  Equation  (210),  produces 


Q  \  2  I  1,2 

—  1  =a  +b  . 


(215) 


Substituting  Equation  (214)  into  Equation  (215)  we  obtain 


Q  \  2 


1  + 


C  sin(6'a: 


Q  -  Ccos{9^) 


=  a 


D 


Q  -  C'cos(6'j,)_ 


(216) 


Solving  Equation  (216)  for  a  yields 


a  = 


Q  {Q  -Ccos{9^)) 
D2 


(217) 


Equation  (217)  is  related  to  the  real  part  of  the  viscoelastic  model  in  terms  of  ex¬ 
perimental  data.  We  may  now  substitute  Equation  (217)  into  Equation  (214)  to 
obtain 


CQ  sm{9x) 


(218) 


Equation  (218)  is  related  to  the  imaginary  part  of  the  viscoelastic  model  in  terms  of 
experimental  data. 

Up  until  this  point,  no  specihc  material  model  has  been  imposed,  i.e.,  A(a;)  and 
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fi{u)  were  related  through  Poisson’s  ratio  to  E{u).  This  relationship  is  stated  in 
Equations  (175)  and  (176).  From  the  discussion  in  Section  3.1,  E{u)  can  take  on 
many  forms,  which  depend  on  the  behavior  of  the  material.  In  order  to  obtain 
meaningful  information  from  Equations  (217)  and  (218),  we  will  choose  to  repre¬ 
sent  the  viscoelastic  model  as  a  general  polynomial  model  as  in  Equation  (52).  Thus, 
let  E{u)  =  p{uu2)  -l-  iq{uu2).  Therefore,  E{1)  =  p{u2)  +  iq{uj2).  We  note  that 
E{0)  =  p(0)  -|-  iq{0).  Recall  from  Section  4.12  that  we  require  E{0)  =  1,  which  leads 
to  p(0)  =  1  and  g(0)  =  0.  Therefore,  equating  a  and  b  to  the  general  viscoelastic 
model  gives 

a  =  l  +  2x*p{oo2),  (219) 


and 


b  =  2x*q{uj2). 


(220) 


Using  Equations  (217)  and  (219)  we  can  solve  for  p{u)2)  to  obtain 


p{uj2)  = 


CQ  cos{e^)  - 
2x*D‘^ 


(221) 


In  addition,  using  Equation  (218)  and  (220),  we  can  solve  for  q{u2)  to  obtain 


?(^^2) 


CQ  sm{9x) 
2x*D‘^ 


(222) 


Equations  (221)  and  (222)  represent  the  solution  to  the  viscoelastic  model  parameters 
in  terms  of  experimental  data.  Keep  in  mind  that  we  require  p(0)  =  1  and  q{0)  =  0. 
We  will  address  these  issues  in  the  following  section  but,  before  we  do  so,  we  must 
hrst  discuss  the  experimental  data. 
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5.4  Further  Analysis  and  Comparison  to  the  Experimental  Data 

Now  that  we  have  the  solution  of  the  model  parameters  related  to  the  experimental 
data,  Equations  (202),  (221),  and  (222),  we  can  use  the  experimental  data  to  generate 
useful  information  about  the  nanoscale  surface  material  properties.  We  shall  analyze 
the  model  and  the  experiment  in  two  pieces.  These  are  the  static  piece,  and  the 
dynamic  piece. 


5.4.1  Static  Piece. 


First,  we  shall  concentrate  on  y*,  which  comes  from  the  static  portion  of  the 
experimental  data.  In  the  experiment,  Qdc  was  not  measured  so,  we  will  have  to 
assume  that  we  know  one  of  the  model  parameters.  In  this  case,  the  material  in  the 
experiment  is  polystyrene.  Polystyrene  has  a  published  value  for  the  Young’s  modulus 
of  about  3  GPa,  and  a  Poisson’s  ratio  of  about  Thus,  the  plane  strain  Young’s 
modulus  for  polystyrene  is  Eq  =  3.375  GPa. 

Under  the  assumption  that  Eq  is  known,  we  will  use  Equation  (201)  to  solve  for 
the  unknown  parameter,  c°,  to  obtain 


c 


0 

a 


3RkL 

^9aEQ 


(223) 


where  we  have  used  the  fact  that  L  =  Next,  substitute  Equation  (223)  into 


X*  =  fo  obtain 


X  = 


3RL 

3 

r^oi 

1 

_ 1 

.  k  _ 

(224) 


This  allows  us  to  estimate  x*  given  the  surface’s  plane  strain  Young’s  modulus,  Eq. 

Knowing  the  effective  cantilever  stiffness  for  the  experiment,  k  =  0.06  N/m,  the 
tip  curvature  radius,  i?  =  30  nm,  and  the  preload  voltage,  L  =  0.25  V,  we  can  use 
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Table  2.  x*  various  preloads. 


L 

X* 

0.25  V 

114.34 

0.5  V 

144.07 

1.0  V 

181.51 

2.0  V 

228.69 

3.0  V 

261.78 

Equation  (224)  to  estimate  x*  =  114.34.  In  the  experiment,  various  preloads  that 
range  in  value  from  .25  V  to  5.4  V  were  used.  Table  2  lists  several  estimated  values 
of  X*- 

We  see  from  Table  2  and  from  Equation  (224)  that  as  the  preload  is  increased, 
X*  also  increases.  Since  the  effective  cantilever  spring  constant  and  the  plane  strain 
Young’s  modulus  of  the  material  remain  constant,  is  directly  proportional  to  x* 
through  c°  =  X*;^-  Thus,  as  we  increase  our  preload  voltage,  the  average  contact 
area  radius  must  increase.  This  is  precisely  what  we  would  expect  to  happen  when 
we  apply  a  greater  DC  loading. 

5.4.2  Dynamic  Piece. 

Now  that  we  have  estimated  the  values  of  x*  for  various  preload  voltages,  we  can 
move  on  to  looking  at  the  viscoelastic  material  model  parameters.  First,  we  must 
examine  the  data  taken  from  the  experiment.  This  data  is  represented  in  terms  of 
the  error  amplitude,  P,  and  phase,  0p,  as  well  as  the  piezo  signal  amplitude,  Q,  and 
phase,  Gq.  Three  sets  of  data  were  collected.  Each  set  of  data  corresponds  to  a 
different  thickness  of  the  polystyrene  surface.  The  surfaces  were  30  nm,  70  nm,  and 
220  nm  thick.  All  of  the  following  results  remain  consistent  over  any  of  the  surface 
thicknesses.  For  this  reason,  we  will  only  present  our  results  for  the  220  nm  thick 
substrate. 
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Figure  7.  The  top  plot  is  the  modulation  frequency  in  Hertz  ver¬ 
sus  the  amplitude  in  Volts  of  the  piezo  drive  signal.  The  bottom 
plot  is  the  modulation  frequency  in  Hertz  versus  the  phase  in  de¬ 
grees  of  the  piezo  drive  signal.  The  experimental  data  was  taken 
with  a  preload  of  0.25  V  on  a  220  nm  thick  polystyrene  surface. 

The  data  was  collected  by  setting  the  AFM  tip  on  the  surface,  and  applying  a  DC 
voltage.  Since  we  are  dealing  with  a  viscoelastic  surface,  there  is  an  initial  creep  before 
‘near’  equilibrium  is  reached.  During  the  viscoelastic  creep,  hve  to  six  data  points 
are  collected.  A  data  point  is  a  measure  of  amplitude  in  volts,  or  phase  in  degrees. 
After  equilibrium  is  reached,  the  sinusoidal  modulation  voltage  bias,  7  =  0.05  V,  is 
applied  at  a  frequency  of  1000  Hz.  The  system  is  allowed  to  reach  steady  state,  then 
another  data  point  is  collected.  After  which,  the  frequency  is  lowered  and  the  process 
is  repeated.  This  is  done  in  a  sweep  from  1000  Hz  to  1  Hz  in  801  equal  increments. 
After  the  last  data  point  is  recorded,  a  new  preload  is  applied,  and  the  frequency 
sweep  is  repeated.  This  is  done  for  10  to  15  preloads,  which  range  from  0.25  V  to  5.4 
V. 

So,  we  must  first  eliminate  the  first  five  to  six  data  points  from  our  data  set 
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Figure  8.  The  top  plot  is  the  modulation  frequency  in  Hertz 
versus  the  amplitude  in  Volts  of  the  error  signal.  The  bottom 
plot  is  the  modulation  frequency  in  Hertz  versus  the  phase  in 
degrees  of  the  error  signal.  The  experimental  data  was  taken 
with  a  preload  of  0.25  V  on  a  220  nm  thick  polystyrene  surface. 

because  this  is  when  the  viscoelastic  creep  is  occurring.  After  that,  the  data  set  is 
divided  into  individual  frequency  sweeps.  This  yields  a  set  of  data  for  each  preload. 
Examples  of  the  data  for  a  220  nm  thick  polystyrene  layer  with  a  preload  of  0.25  V 
are  given  in  Figures  7  and  8. 

Examining  the  amplitude  and  phase  of  the  piezo  drive  signal  in  Figure  7,  we  see 
that  as  the  driving  frequency  increases,  there  is  a  slight  initial  increase  in  amplitude. 
This  is  followed  by  an  almost  linear  decrease  in  amplitude  starting  from  about  200 
Hz.  The  quick  drop  offs  at  near  zero  frequency  and  at  near  1000  Hz  are  the  result 
of  the  changeover  between  the  differing  frequency  data  sweeps.  The  phase  of  the 
piezo  is  seen  to  steadily  increase  as  the  tapping  frequency  is  increased.  The  small 
jump  in  phase  at  around  950  Hz  is  an  as  yet  unexplained  anomaly.  Turning  to  the 
error  data  in  Figure  8,  we  see  that  at  near  zero  frequency,  there  is  almost  no  error  in 
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Figure  9.  These  plots  are  calculated  results  from  our  model. 

The  top  plot  is  the  modulation  frequency  in  Hertz  versus  the 
amplitude  in  Volts  of  the  AFM  tip  position.  The  bottom  plot  is 
the  modulation  frequency  in  Hertz  versus  the  phase  in  radians  of 
the  AFM  tip  position. 

the  amplitude.  This  means  that  at  low  frequency,  the  AFM  has  the  same  amplitude 
as  the  driving  amplitude.  As  the  frequency  is  increased,  there  is  a  steady  increase 
in  the  error  amplitude.  The  phase  near  zero  frequency  is  about  90  degrees.  As  the 
driving  frequency  is  increased,  the  phase  steadily  decreases.  Again,  we  note  the  as 
yet  unexplained  anomalous  jump  around  950  Hz. 

Now,  we  are  going  to  calculate  the  dynamic  cantilever  tip  position  amplitude, 
(7,  from  Equation  (194).  We  will  also  calculate  the  phase,  9c,  of  the  cantilever  tip 
position  from  Equation  (195).  Plots  of  the  tip  position  amplitude  and  phase  are 
shown  in  Figure  9. 

From  Figure  9  we  can  see  that  as  we  approach  zero  frequency,  the  AFM  tip  moves 
at  the  driving  amplitude  and  phase.  Exactly  as  we  would  expect  since,  the  material 
response  time  will  be  faster  than  the  driving  frequency.  This  means  that  the  surface 
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Figure  10.  This  plot  is  calculated  results  from  our  model.  This 
plot  is  the  modulation  frequency  in  Hertz  versus  the  amplitude 
in  Volts  of  the  AFM  tip  position  for  various  preload  voltages. 

responds  fast  enough,  and  the  cantilever  moves  in  phase  and  at  the  same  amplitude 
as  our  modulation  voltage.  As  the  frequency  increases,  we  see  the  viscoelasticity  of 
the  surface  lowers  the  amplitude  of  the  tip  movement.  We  also  see  the  phase  shift 
increases  in  magnitude  as  we  increase  in  frequency. 

We  will  now  present  more  results  for  the  calculated  AFM  tip  position  amplitude, 
C,  and  its  phase,  9c,  based  on  the  data  for  various  preloads  at  220nm.  The  calculated 
AFM  tip  position  amplitude  for  various  preloads  is  shown  in  Figure  10.  At  lower 
frequencies,  the  cantilever  tip  position  is  larger  when  there  is  a  higher  preload.  As  we 
increase  in  frequency,  the  AFM  tip  position  amplitude  begins  to  decrease  regardless 
of  the  preload.  Also,  the  slope  of  the  amplitude  curve  is  greater  for  higher  preloads. 

Looking  at  the  phase  from  the  cantilever  tip  movement,  in  Figure  11,  we  see  that 
at  low  driving  frequencies  the  AFM  tip  is  in  phase  with  the  surface.  As  the  driving 
frequency  is  increased,  the  phase  shift  decreases.  The  larger  the  preload  voltage,  the 
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Figure  11.  This  plot  is  the  modulation  frequency  in  Hertz  versus 
the  phase  in  radians  of  the  AFM  tip  position  for  various  preload 
voltages. 

greater  the  rate  of  decrease  in  the  phase  shift.  At  just  over  600  Hz  for  a  3  volt  preload, 
we  see  a  jump  in  the  phase  and  the  amplitude  of  the  tip  position.  This  is  because  the 
tip  has  such  a  high  preload  and  is  180  degrees  out  of  phase  with  the  driving  voltage. 
We  believe  the  tip  actually  left  the  surface,  which  is  why  there  is  a  jump  in  both  the 
amplitude  and  phase.  Observe  that  the  unexplained  anomalous  phase  jump  around 
950  Hz  presents  itself  in  the  AFM  phase  data  as  well. 

Next,  we  examine  the  viscoelastic  model  parameters.  Since  the  piezo  delay  time 
was  not  measured,  we  can  only  analyze  the  viscoelastic  parameters  from  a  purely 
mathematical  standpoint.  First,  recall  that  we  need  to  satisfy  p{0)  =  1  and  g(0)  =  0 
in  the  viscoelastic  model.  Equations  (221)  and  (222).  Note  from  Figures  7  and  9  that 
both  the  piezo  amplitude,  Q,  and  AFM  amplitude,  C,  are  continuous  with  respect  to 
frequency.  Also,  approaching  zero  tapping  frequency  C,Q  >  0.  Therefore,  as  long  as 
11(0)  7^  0,  then  Equation  (222)  requires  03,(0)  =  0  in  order  to  satisfy  g(0)  =  0. 
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To  demonstrate  that  -D(O)  7^  0,  suppose  not.  For  Z)^(O)  =  0,  using  Equation  (208) 


requires  cos(6'a:(0))  =  '^2C(o|q(o)°^  •  Thus, 


Since  both  C(0)  and  <5(0)  are  real,  then  Equation  (225)  can  only  be  satisfied  if  6^3, (0) 
is  complex.  However,  0^  =  Oc  —  Oq  +  t*  is  a  real  number.  Therefore,  £>^(0)  7^  0. 

Now,  using  the  requirements  that  p(0)  =  1  and  0^(0)  =  0,  along  with  Equa¬ 
tion  (221)  we  obtain 


C'(O) 


(226) 


2(Q(0)-C'(0))- 


This  allows  us  to  estimate  x*  without  having  to  know  Qdc  as  was  necessary  in 
Equation  (202).  The  only  requirement  is  that  t*  approaches  6q  —  6c  near  zero  tapping 
frequency.  Since  the  viscoelastic  material  has  a  faster  wave  speed  than  the  driving 
frequency  near  zero  (see  Figure  9  and  the  argument  given  earlier  in  this  section),  then 
we  also  have  a  way  of  estimating  the  piezo  delay  time  near  zero  tapping  frequency. 

Although  the  preceding  analysis  means  that  the  mathematical  model  is  consistent 
with  the  experiment,  without  the  actual  measurements  of  the  DC  piezo  loading,  Qdc-, 
and  the  piezo  delay  time,  r^,  we  cannot  truly  give  any  more  material  properties.  Since 
Tp  is  a  function  of  the  driving  frequency,  then  knowing  an  estimate  of  the  delay  time 
near  zero  frequency  will  not  help  to  calculate  the  frequency  dependent,  viscoelastic 
material  properties.  In  addition,  the  variability  in  how  the  data  is  collected  makes  it 
very  difficult,  and  possibly  inaccurate,  to  estimate  x*  based  on  Equation  (226).  That 
being  said,  an  estimate  was  attempted  using  a  linear  fit,  and  for  a  preload  of  0.25  V, 
y*  50.  Using  Equation  (224)  this  means  that  Eq  1  GPa,  which  is  of  the  same 
order  of  magnitude  as  our  estimate  based  on  published  values. 


In  summary,  this  chapter  sets  up  an  error  model  that  allows  us  to  calculate  the 


penetration  depth  of  the  AFM  tip  into  the  material  based  on  the  experimental  data. 
We  then  equated  this  solution  to  our  analytical  solution  for  the  penetration  depth 
from  the  last  chapter.  This  allowed  the  production  of  an  analytical  solution  of  the 
viscoelastic  material  parameters  in  terms  of  experimental  data.  We  then  presented 
our  calculated  data  and  made  observations  based  on  our  results.  The  next  chapter 
will  present  our  conclusions  and  other  possible  future  extensions  and  applications  of 
this  analytical  AFM  model. 
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VI.  Conclusions 


Modeling  the  AFM  dynamic  nanoindentation  experiment  that  was  described  in 
Section  1.2  provides  us  with  a  wealth  of  information.  The  model  of  the  experimental 
setup  was  divided  into  two  models  that  are  coupled  by  boundary  conditions.  The  two 
models  created  were  an  AFM  model  and  a  material  model. 

For  the  AFM  model,  we  initially  assume  we  can  model  the  AFM  probe  as  a 
cantilevered  beam.  This  leads  us  to  a  one-dimensional  spring-mass  representation  of 
the  AFM.  The  forcing  functions  used  in  the  spring-mass  system  are  a  static  force,  a 
sinusoidal  force,  and  the  force  generated  by  the  surface  stresses  over  the  area  of  the 
AFM  tip.  The  surface  stress  forces  are  used  to  couple  the  AFM  model  to  the  surface. 
The  static  and  sinusoidal  forces  are  generic  and  in  the  analytical  model  solution 
they  represent  the  piezo  forcing  of  the  system.  It  should  be  noted  that  AFM’s  and 
nanoindentation  systems  can  both  be  modeled  as  spring-mass  systems,  and  with  very 
little  modihcation,  the  analytical  model  produced  from  this  research  can  be  modihed 
to  incorporate  either  system.  This  shows  the  possibility  of  broad  application  of  our 
analytical  model. 

The  material  model  is  based  on  viscoelastic  material  behavior.  The  model  assumes 
an  axisymmetric  indenter  which  indents  a  semi-inhnite  half-space.  The  indentation 
is  assumed  to  be  a  static  displacement  superposed  with  a  much  smaller  dynamic 
displacement.  Specihc  boundary  conditions  are  chosen  to  represent  an  AFM  dynamic 
nanoindentation  experiment  which,  as  it  turns  out,  are  the  same  conditions  as  a  (non- 
AFM)  nanoindentation  experiment.  The  analytical  solution  of  the  material  model  for 
the  penetration  depth  of  the  indenter  contains  a  static  and  a  dynamic  piece.  The  static 
piece  corresponds  directly  to  the  viscoelastic  extension  of  Sneddon’s  solution  [52, 
84],  which  lends  us  to  believe  the  dynamic  piece  is  correct.  Current  solutions  only 
modify  the  viscoelastic  extension  of  Sneddon’s  solution  in  light  of  small  dynamic 


displacements.  This  violates  certain  assumptions  in  the  problems  formulation,  and  is 
mathematically  incorrect. 

The  coupled  AFM  and  material  models  are  then  analytically  solved  for  the  pen¬ 
etration  depth  of  the  AFM  tip  into  a  viscoelastic  surface.  The  viscoelastic  model  is 
assumed  to  be  represented  by  a  polynomial  in  frequency,  and  can  be  chosen  after  the 
experimental  data  is  analyzed.  Current  solutions  assume  a  viscoelastic  model  before¬ 
hand,  and  are  then  pigeonholed  by  the  formulation  of  the  model.  It  should  be  noted 
that  the  solution  for  the  penetration  depth  is  under  the  assumption  of  forcing  at 
low  frequencies.  The  solution  can  easily  be  modified  to  incorporate  forcing  at  higher 
frequencies. 

A  third  model  we  called  the  error  model  was  then  created  to  relate  the  raw  exper¬ 
imental  data  to  the  penetration  depth  of  the  AFM  tip.  The  analytical  solution  for 
the  penetration  depth  of  the  coupled  AFM  and  material  model  was  then  equated  to 
the  error  model  solution.  This  allowed  us  to  produce  a  simple  analytical  model  that 
relates  raw  experimental  data  to  viscoelastic  material  parameters.  The  manipulated 
data  can  then  be  fit  with  a  polynomial  model,  and  the  coefficients  of  the  fit  can  be 
interpreted. 

As  an  illustration,  viscoelastic  properties  for  a  polystyrene  thin  him  on  a  silicon 
substrate  are  calculated  and  presented  using  the  analytical  solution.  Because  of  lack  of 
experimental  measurements,  not  all  calculations  were  possible.  It  is  our  assertion  that 
given  all  the  necessary  experimental  data,  our  simple  analytical  model  will  provide 
viscoelastic  properties  of  near-surface  materials  from  AFM  dynamic  nanoindentation 
experiments.  The  novelty  of  our  analytical  model  is  that  it  admits  different  stress- 
strain  materials  descriptions  than  current  models. 

In  addition  to  the  extensions  in  applications  of  this  model  to  broader  experiments, 
as  described  above,  the  model  may  be  extended  in  other  ways.  Possible  extensions 


include  the  incorporation  of  adhesion  to  the  material  model,  though  this  will  likely 
produce  integrals  with  only  numerical  solutions.  Other  possibilities  include  the  cal¬ 
culation  of  the  next,  smaller  term  in  the  linear  approximation  of  the  material  model. 
The  next  term  will  contain  convolutions  of  functions  and  in  light  of  the  currently 
modeled  experiment  was  considered  negligibly  small. 

This  research  in  nanoscale  viscoelastic  material  properties  allows  for  the  advance¬ 
ment  of  study  in  nanomaterials  and  has  applications  in  many  areas.  Of  particular 
interest  to  the  United  States  Air  Force  and  the  Department  of  Defense  are,  high- 
altitude  long-endurance  ISR  airships,  prompt  theater- range  ISR/strike  systems,  direct 
forward  air  delivery  and  resupply,  energy-efficient  partially  buoyant  cargo  airlifters, 
fuel-efficient  hybrid  wing-body  aircraft,  and  hyperprecision  low-collateral  damage  mu¬ 
nitions  [17].  This  research  will  help  further  all  helds  that  are  associated  with  nanome¬ 
chanical  properties. 


90 


Appendix  A.  Cylindrical  Coordinates 


A 


Figure  12.  The  cylindrical  coordinate  system. 


The  cylindrical  coordinate  system  is  a  generalization  of  polar  coordinates  to  three 
dimensions.  The  representation  of  a  point  in  cylindrical  coordinates  is  {r,6,z),  as 
can  be  seen  in  Figure  12.  Simple  formulas  exist  to  convert  from  the  rectangular 
coordinates  {x,y,z)  to  the  cylindrical  coordinates  {r,9,z).  These  formulas  are 

r  =  ^x^  +  y^,  (227) 

6  =  arctan  (—')  ,  (228) 

\  X  / 

and 

z  =  z,  (229) 

where  r  G  [0,oo),  9  G  [0,27r),  and  z  G  (—00,00).  The  inverse  tangent  function  must 
be  dehned  in  order  to  take  the  correct  quadrant  of  {x,  y)  into  account. 
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In  order  to  convert  from  the  cylindrical  coordinates  {r,6,z)  to  the  rectangular 


coordinates  {x,y,z),  we  use  the  formulas 


and 


X  =  r  cos  {6) , 

(230) 

1/  =  r  sin  {6) , 

(231) 

z  =  z.  (232) 

The  unit  vectors  in  cylindrical  coordinates  are  f ,  6,  and  z.  They  dehne  the  direc¬ 
tions  of  the  coordinate  terms,  and  can  be  seen  in  Figure  12.  We  dehne  the  gradient 
as 


_  ^  d  -Id  ^  d 


(233) 


The  divergence  is  dehned  as 


V  •  u 


r  dr 


{rur)  + 


1  due 
r  do 


+ 


du^ 
dz  ' 


(234) 


The  curl  is  dehned  by 


V  X  M  = 


duB\  f  dUr 


dr  J  r 


{rue) 


(235) 
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Appendix  B.  Dual  Integral  Equations 


This  appendix  shows  a  brief  overview,  without  derivation,  of  the  results  of  Man- 
dal  [60].  For  a  complete  derivation,  please  refer  to  that  paper. 

The  general  form  of  dual  integral  equations  involving  Bessel  functions  of  the  first 
kind  is 

Ju{pf3)  [1  +  w{l3)]  $(/3)d/3  =  /(p),  0  <  p  <  1, 

Jf,{pl3)^{l3)dl3  =  g{p),  p  >  1,  (236) 

where  a*  and  (3*  are  known  constants,  v  and  p  are  the  known  orders  of  the  Bessel 
functions,  tc(/9)  is  an  arbitrary  weight  function,  /(p)  and  g{p)  are  known  functions 
valid  over  their  particular  region,  and  *F(/3)  is  the  unknown  function  to  be  determined. 

Through  the  use  of  Sonine’s  integrals  and  Hankel  inversion,  Mandal  reduced  Equa¬ 
tion  (236)  to  a  Fredholm  integral  equation  of  the  second  kind.  This  equation  in  $(p) 
is 


$(p)  =  p 


_  7-s+l 


••I  noo 

rFi{r)Jx+s{rp)dr  +  /  rGi{r)Jx+s{rp)dr 
L-'o  Ji 

-7+s  /  \'^J\+s+l{u)Jx+s{,p')  ~  pJ\+s+l{,P)J\+s{u)  , 

u  '  t(;(M)<I)(M) - - - r - du 

—  p^ 


,  (237) 


where  A  and  7  are  defined  by 


7  = 


2 

p  —  p 


+  13*  +  a\ 


(238a) 

(238b) 


and  s  is  a  nonnegative  integer  we  choose  so  that  A  +  s  >  —1.  Fi(r)  and  Gi(r)  are 
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functions  that  are  defined  by 


Fi(r) 

Gi(r) 


(-1)V+" 


(239a) 

(239b) 


Here,  the  functions  F{r)  and  G{r)  are  dehned  as 

('—1 /I  ^  \  roo 

GW  =  (;^)  I  ^-'■^Yp^-rTMdp,  (240b) 

where  r(2;)  is  the  Gamma  function  and 


e  -  a*  -  1  +  £,  (241a) 

T]  =  y-  +  a*  —  /?*  —  1  +  m.  (241b) 

The  terms  i  and  m  are  nonnegative  integers  that  we  choose  so  that  ^  >  —  1  and 
rj  >  —1.  The  requirements  that  certain  terms  are  greater  than  —1  are  because  of  the 
Hankel  inversions  used  to  solve  Equation  (236). 
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Appendix  C.  Integral  Equation  Neumann  Series  [41] 

An  integral  equation  Neumann  Series  is  used  to  solve  Fredholm  integral  equations 
of  the  second  kind.  A  Fredholm  integral  equation  of  the  second  kind  takes  the  form 

(j){x)  =  f{x)  +  ^  f  K{x,t)(j){t)dt,  (242) 

J  a 

where  a,  6  G  M  and  A  G  M  —  {0}.  Dehne  the  operator  K  by 

K(j){x)  =  f  K{x,t)(p{t)dt.  (243) 

J  a 

Therefore,  we  can  write  Equation  (242)  as 

Tct>  =  0,  (244) 

where 

T(f)  =  f  +  XK(j).  (245) 

If  /  is  in  a  Hilbert  space  H,  and  K  is  a  bounded  linear  operator  with  the  property 


||iF0i-iF02||  <M||0i-02||,  (246) 

then  Equation  (244)  has  a  unique  solution  for  all  /  provided  that  |A|M  <  1.  In 
addition,  if  /  is  in  L2[a,fe],  then 

0(x)  =  lim  TVo(^),  (247) 

n^oo 

where  fo{x)  G  L2[a,h]  is  an  arbitrary  initial  function.  So, 

T/o  =  /  +  AiF/o,  (248) 
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TVo  =  T[f  +  XKfo]  =  f  +  \Kf  +  X^K^fo, 


(249) 


and 

T"/o  =  /  +  XKf  +  X^K^f  +  ■  ■  ■  +  X^-^K^-^f  +  X^K^fo.  (250) 


Therefore, 


0  =  /  +  XKf  +  X^K^f  +  ■  ■  ■  +  +  ■  ■  ■  .  (251) 


We  define  the  integral  operator  as 


^”/(^)  =  /  Kn{x,y)f{y)dy,  (252) 

J  a 


where  Kn{x,y)  is  dehned  recursively  by 


Kn{x,y)=  f  K{x,z)Kn-i{z,y)dz,  n  =  2,3,...,  (253) 

J  a 


and 


Ki{x,y)  =  K{x,y). 


(254) 
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Appendix  D.  Some  Integrals 


This  appendix  includes  some  important  integrals  and  their  derivations. 


PMPp)jo(,P)dP  = 


■  /3)]jo{u)du  =  Tijoi/d). 

(255) 

■  [^)]j2iu)du  =  7lj2{l3). 

(256) 

{1-pTK  0<p<1, 

(257) 

0, 


1  <p. 


I  (1  -  2  (2  -  3p2) ,  0  <  p  <  1, 

/  l3Jo{l3p)j2mi3  = 

'0  10,  1  <  p. 


(258) 


MP)jo{P)d/3  =  1. 
Ji{P)j2{P)dp  =  0. 


(259) 

(260) 


The  derivation  for  Equation  (255)  is  as  follows 


[jo{u  +  /3)+  jo{u  -  P)]  io{u)du 


sin(M  +  /3)  ^  sin(M  —  (3) 


u  +  (3 


u  —  (3 


sin(M) 


du 


u 


sin(M)  cos(/3)  +  sin(/3)  cos(m)  ^  sin(M)  cos(/3)  —  sin(/3)  cos(m) 


=  2 


u  +  (3 

u  sin(M)  cos(/3)  —  (3  sin(/3)  cos(m) 
—  (3“^ 


smfMj 

u 


u  —  (3 
du 


/  sm  (m)  ,  n  r  /  sm(M)  cos(m)  , 


=  2cos(/3) 


sin^iu) 


du  —  (3  sin(/?) 


sin(2M) 

/q  u{u^  —  (3‘^) 


du 


sin(M) 


du 


u 


(261) 
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Now,  from  Gradshteyn  and  Ryzhik  [27]  page  432  we  have 


and  from  page  425  we  have 

I  11  -  “^(2/5)1  =  ^  (263) 


So, 


bo(M  +  /3)  +  Jo(m  -  /5)]  3o{u)du 


=  ^sin(/?)cos^(/3)  +  ^sin3(/?) 


(264) 


=  ^sin(/3)  [cos^(/3)  +  sin" 
P 


The  derivation  for  Equation  (256)  is  as  follows 


sin(/3)  .  ,  . 

=  =  '^doW)- 


[jo{u  +  (3)+  jo{u  -  /3)]  j2{u)du 
=  2 


poo 

usm{u)  cos(/3) 

—  (3  sin(/3)  cos(m) 

^sin(M)  ^cos(m)  sin(M) 

L 

u^ 

-(3^  \ 

3  „  ^2 

U'^  u 

du 


=  6  cos(/3) 


sin^(M)  f°°  sm{u)  cos{u) 

[J,  Jo  n  («2  -  (32) 


du 


6/3  sin(/?) 


•  sin  M  cos  M 


du  + 


cos^(-u) 


Jo  {u^  —  (3"^)  Jo  ~ 


-du 


(265) 


Now,  using  partial  fraction  decomposition,  we  may  write 

I  , 

-du  + 


sin^fw) 


°°  —  sin^fw) 


°°  sin^fw) 


0 


du. 


(266) 
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From  Gradshteyn  and  Ryzhik  [27]  page  431  we  have 


p -sin2(n)^  _  -TT 

io 


(267) 


and  from  page  432  we  have 


sin^(M)  , 

_ (111 

/32  (m2  _ 


TT 


sin(2/3). 


(268) 


Now,  we  may  write 


sin(M)  cos(m)  /■°° 

/,  =  ' 


,  sin(2M)  ,  TT  ,  / 


(269) 


where  we  have  used  Gradshteyn  and  Ryzhik  [27]  page  425.  So, 


sin2(M)  ,  p 
-du  — 


sin(M)  cos(m) 
/g  {u^  —  (5“^)  Jq  u  {u^  —  (5“^) 


du 


-^[l  +  cos(2^)]  +  ^sin(2^) 
[cos2(^)]  +  ^sin(/3)cos(/3). 


(270) 


Now,  by  partial  fraction  decomposition  we  have 


sm  M  cos  m 


/g  {u^  —  /3^) 


du 


smi^Mj  COSl^Mj 

2/1,3 


(3‘^u 


du  + 


sm  M  cos  M 


/3^-u 


du 


MSm  M  cos  M 


/g  (m2  -  ^2) 


dw.  (271) 


Now,  from  Gradshteyn  and  Ryzhik  [27]  page  431  we  have 


Sm(Mj  COS(m)  TT 

*  =  W 


(272) 
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and  from  page  424  we  have 


Msin(M)  cos(m) 
(^2  _  ^2) 


du 


(‘°°  usm{2u) 

io  WW^) 


71 


cos(2/3). 


(273) 


Thus, 


1'°°  —  sin(M)  cos(m) 

io 


SinlMl  COSlMl  ,  TT  ,  / 

-^j^du Mm 


sin(2M) 

2/32^3 


(274) 


Now,  we  may  write 


cos^(m) 

m2  {u^-i32f^ 


cos^(m) 


du  + 


cos^(m) 

^2 


(275) 


From  Gradshteyn  and  Ryzhik  [27]  page  461  we  have 


cos^(m) 

P2 


TT 

4^ 


sm(2;5) 


^sm(/3)cos(/3). 


(276) 


Thus, 


cos^(m) 

m2 


I  +  cos(2m) 

I 


TT 


sin(/?)  cos(/3). 


(277) 


Therefore, 


-sm(M)cos(M)  P  cos^(m) 

m3  (m2  -  /32)  +  („2  _  ^2) 

1  /■°°  sin(2M)  —  u  —  u  cos(2m) 

= 


^  sin(/3)  cos(/3)  +  ^  sm2(/3). 

(278) 
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Now,  we  may  write 


sin(2M)  —  u  —  u  cos(2m) 


du 


1 

2^2  _ 
1 


sin(2M)  —  2mcos(2m)  ,  f°°  mcos(2m)  —  u  , 

— ^ ^ 5 ^ — -du  +  /  du 


'0 
poo 


Jo 

sin(2M)  —  2m  cos(2u)  /“°°  sin^fw) 

- - - du  —  2  - - — du 


uo 


'0 


u^ 


(279) 


From  Gradshteyn  and  Ryzhik  [27]  page  447  we  have 


sm(2M)  —  2mcos(2m) 


du  =  TT, 


(280) 


and  from  page  431  we  have 


r  sm\u)  , 

-2  /  - ^du  =  -71. 


>0 


u^ 


(281) 


Therefore, 


sin(2M)  —  u  —  u  cos(2m) 


2^^ 


du  =  0. 


(282) 


Hence, 


sm  M  cos  M 


cos^(m) 


Iq  U^  (m2  — 


*+  L 


(283) 
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Now,  combining  Equations  (265),  (270),  and  (283)  we  get 


[jo{u  +  /3)  +  io{u  -  (3)]  j2{u)du 


=  6  cos(/3) 

+  6/3  sin(/3) 
sin(/3) 


TT 


2^2  ^  sin(/3)  cos(/3) 


TT 


TT 


2/33 


=  TT 


=  TT 


sin(/3)  cos(/3)  +  ^sin  (/3) 
cos(/3) 


-  7rjo(/5) 


/33 


(cos2(/3)  +  sin2(/3))  -  3^^  (cos2(/3)  +  sin2(/3)) 


gSiii(/3)  _  g  Cos(^)  sin(/3) 


/33 


/32 


/32 

=  7rj2(/3). 


sin(/3) 


(284) 


The  derivation  for  Equation  (257)  is  as  follows 


’  sinf 

/3Jo(/3p)jo(/3)d/3  =  /  (3Jom^^d/3 

Jo  P 

poo 

=  /  Jo(/3p)sin(/3)d^. 


(285) 


Now,  from  Gradshteyn  and  Ryzhik  [27]  page  718  we  have 


I  (1  -  p2)  %  0  <  p  <  1, 

Jo(/3p)  sin(/3)ci/3  = 

'0  10,  1  <  p. 


(286) 


Thus, 


I  (1  -  p^)  G  0  <  p  <  1, 

pjomump  = 

'0  10,  1  <  p. 


(287) 


The  derivation  for  Equation  (258)  is  as  follows 


/3^o(/3p)j2(/3)d/3  =  /  /3Jo(/3p) 


d/3,  (288) 
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where  we  have  used  the  recurrence  relationship  from  Abramowitz  and  Stegun  [1]  of 


jn-l(z)  +jn+liz)  =  {2n+l)z 


(289) 


with  n  =  1.  Now,  using  the  fact  that  jn{x)  =  we  can  write  Equa¬ 

tion  (288)  as 


f3JoWp)j2{f3)df3  =  ^  /  13  ^/Vo(^p)J3/2(/3)d^  -  /  f3Jo{f3p)jo{l3)df3. 


From  Gradshteyn  and  Ryzhik  [27]  page  683  we  have 


(290) 


/3-V2jo(/3p)  J3/2(/5)d/3  = 


0, 


^(l-p2)^  0<p<l, 

1  <  p. 


(291) 


Finally,  combining  Equations  (257)  and  (290)-(291)  yields 


■  (2-3p2)(l-p^)-%  0<p<l, 

/  l3Jo{l3p)]2{l3)d(3  = 

'0  10,  1  <  p. 


(292) 


The  derivation  of  Equation  (259)  is  as  follows 


Mfi)Mfi)dD  =  r 


(293) 


Now,  from  Gradshteyn  and  Ryzhik  [27]  page  727  we  have 


jdPf-P^dp  =  i. 


(294) 


So, 


Ji{(3)jo{/3)df3  =  1. 


(295) 


103 


The  derivation  of  Equation  (260)  is  as  follows 


(296) 


where  we  have  used  the  fact  that  in{,x)  =  ^/^Jn+i/2{x).  Now,  from  Gradshteyn  and 
Ryzhik  [27]  page  683  we  have 


=  0. 


(297) 


Thus, 


Ji{P)j2{P)dP  =  0. 


(298) 
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